For Reference 


NOT TO BE TAKEN FROM THIS ROOM 


Gx AIBRIS 
ANINTRSUTARIS 
aE RUAKASIS 


The University of Alberta — 


| 
Printing Department 
Edmonton, Alberta | 


_ Digitized by the Internet Archive 
in 2023 with funding from 
University of Alberta Library 


https://archive.org/details/Alpaslan19 74 


THE UNIVERSITY OF ALBERTA 


TRANSFER FUNCTION ANALYSIS 


by 


TUMER ALPASLAN 


A THESIS 
SUBMITTED TO THE FACULTY OF GRADUATE STUDIES AND RESEARCH 
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE 


OF DOCTOR OF PHILOSOPHY 


DEPARTMENT OF PHYSICS 


_ EDMONTON, ALBERTA 


SPRING 1974 


ABSTRACT 


The spectral ratio method is used to investigate 
the transmission characteristics of the crust under a 
group of seismic stations in Alberta. Various model 
crustal spectral ratios using a Haskell-Thomson matrix 
formulation have been compared with the experimental 
ones. 

The implications of a simple form of seismic 
anisotropy in plane parallel layered media have been 
investigated. In the low frequency band appropriate 
to the P coda of most teleseismic events, transverse 
anisotropy does not have a significant effect on 
spectral ratios. A ray theory approach is used to 
evaluate the response of a multidipping layered media. 
A computerized method is used to obtain the number of 
terms necessary in a partial ray eRAA Se HOHE approximate 
the wave solution for a flat layered case. Model studies 
indicate that a flat layered assumption is valid for 
dip angles less than 10°. However, low velocity 
layers dipping as much as 10° introduce significant 
shift in peak positions. Complex velocities are in- 
troduced into a Haskell-Thomson matrix formulation to 
study the attenuation characteristics of the crust. 
Comparison of experimental spectral ratios with the 
attenuating spectral ratio model favors an average Qp 


value of 500 in the crust. 
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Analysis of seismic station recordings near 
Edmonton, Pine Lake, and Rocky Mountain House suggests 
_the presence of an approximately 35 km thick crust. 
The station near Trochu and Delburne indicates a 
very heterogeneous structure. Azimuthal deviation 
meastinenente and record character suggest that the 
anomalous results are due to steeply dipping or more 
complicated structure in the upper crust. 

Synthetic seismograms for attenuating isotropic 
or anisotropic layered media have been obtained for an 
input plane P or S wave by taking the inverse Fourier 
transform of the transfer function calculated by a 
Haskell-Thomson matrix formulation. Model seismograms 
have also been calculated for dipping layered media 
using asymptotic ray theory. When an S wave is inci- 
dent at the bottom of the layered media, conversions 
at various interfaces cause a series of precursors to 
S wave. On the vertical component of displacement, 
the main "S" arrival may be 1-6 seconds later than the 
precursors for a continental crust. At some epicentral 
distances a sedimentary layer may generate very strong 
precursors up to 2 seconds before the pure shear 
arrival. This result appears to explain many inconsis- 


tent interpretations in the study of shear waves. 
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CHAPTER I 


INTRODUCTION 


Xe Object of this work 


In seismology a study of the transmission of 
elastic waves has been one of the challenging research 
areas. Tie Beano of Geophysics offers a direct or 
indirect measurement of parameters of the Earth's 
interior. Seismological measurements have shown that 
the interior of the Earth can be divided into several 
distinct physical and chemical shells. The outer part 
of the Earth which is about 50 km thick under continents 
and 5 km under the oceans has Been named the crust. 
This region is separated from the upper mantle by a 
well recognized discontinuity called the "Moho". Rapid 
variations in the properties of the crust are also well 
known. Mountains anid their roots, deep oceans and other 
complex tectonic features create difficulties when 
interpreting elastic waves that have passed through 
the deeper parts of the Earth. 

It is useful sometimes to express the recorded 
teleseismic signal as a convolution of a pulse with a 
transfer function representing the transmission charac- 
teristics of the crust near the recording seismic 
station. There are a number of theoretical techniques 


for calculating the transfer function of layered 
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media. These methods have been extended to yield 
complete wave solutionsfor an incident plane wave 

on a plane parallel layered media in which seismic 
anisotropy and attenuation may exist. In this thesis 
the problem of a dipping layered medium has also been 
investigated independently. The models used in this 
work are assumed to be homogeneous and the layers 

are separated by first order discontinuities. A 
comparison between theory and experiment has been 
attempted. In this way it is hoped that the trans- 
mission characteristics of any crustal section may be 
determined. A number of synthetic seismograms and 
spectral models provided in this work should help 
seismologiststo understand the wave modifications that 
are produced during transmission of a wave through a 


layered medium. 


2. Previous work 


Observations on the effect of the crust on wave 
transmission has proceeded since the early days of 
seismology. Among these Imamura (1929) studied the 
vibration of the crust using seismic waves from earth- 
quakes; Nasu (1931) investigated the relation between 
the period of seismic waves and the type of soil; 
Gutenberg (1957) also studied earthquake damage and 


its relation to the type of soil; Nuttli and Whitmore 
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(1961, 1962) investigated the polarization angle 
between P and S waves and suggested frequency depen- 
dent polarization angle. 

A number of authors attempted to calculate the 
theoretical response of crustal structure. The crust 
is assumed to consist of a number of layers. In the 
early days, great simplifications were made to obtain 
a solution [Sezawa and Kanai’ (1932, 1937)]. 

Thomson (1950) and Haskell (1953) were the first to 
formulate the frequency dependent polarization angle 
characteristics of a layered medium. They expressed the 
response of a plane parallel isotropic layered medium 
to an incoming plane wave at any angle of incidence 
using products of four by four matrices whose elements 
are’ a~ function” of* parameters: of™each’Layer” and*the boundary 
conditions. Dorman (1962) formulated the solution at 

a boundary between a liquid and a solid. Harkrider and 
Anderson (1962) obtained the response for a transversely 
isotropic medium. 

Fast digital computers have made it practical 
to use the Haskell-Thomson matrix method for modelling 
a crustal section with any number of layers. Phinney 
(1964) compared the ratio of long-period spectra of 
(eeTGb aa and horizontal seismograms with a number of 


models calculated using the Haskell-Thomson method. 
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The introduction of a spectral ratio technique was a 
significant breakthrough because it eliminated most 

of the source function effects. Hannon (1964) first 
attempted the calculation of long period synthetic 
seismograms by taking the inverse Fourier transform 

of the transfer function. Fernandez (1967) used 
spectral ratios for long-period data giving model 
curves to determine crustal structure. Leblanc (1967) 
discussed the importance of a window for the high 
frequency portion of the spectrum using une spectral 
ratio method. McCamy (1967) attempted to apply transfer 
function analysis to determine crustal structure using 
everson sources. Ellis and Basham (1968) applied the 
spectral ratio method in an attempt to resolve crustal 
structure in central Alberta. Hasegawa (1970,1971) 
used the spectral ratio method to model the structure 
at Yellowknife. Ibrahim (1969) made use of the spectral 
behavior of SH waves to determine crustal thicknesses. 
Kurita (1969) used long-period waves in an attempt to 
find crustal and mantle structure. Ishii and Ellis 
(1970) formulated the response of single dipping layer 
over a half space and Rogers and Kisslinger (1972) 
applied this ray theory to earthquake data. Jensen 
(1970) applied linear filter theory to calculate the 
response of a layered medium in the time domain to 


reduce the inefficiency of the Haskell-Thomson method. 
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Sprenke (1972) studied the transmission characteristics 
of the Alberta crust using the Haskell-Thomson method. 
In this study the Haskell-Thomson matrix method 
and a ray theory approach were applied to earthquake 
data recorded by University of Alberta in 1970. Impli- 
cations of an anistropic or dipping crustal model are 
discussed. Comparisons are made between the completed 
wave solution and a partial ray expansion and criteria fora 
satisfactory solution have been established both in the 
frequency and the time domain. High frequency theore- 
tical seismograms were obtained by Fourier inversion of 


Haskell-Thomson matrices using a synthetic pulse. 
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CHAPTER II 


TRANSFER FUNCTION STUDIES 


I. Crustal Investigation using the Spectral Ratio Method 


1. Crustal transfer function 


In the time domain, according to filter theory, 
the transmission of a plane elastic wave through a 
mediumcan be represented as a convolution of a pulse 
with a transfer function representing the elastic 
response of the media. In the frequency domain it can be 


written as 
R(w) = P(w) -H(w) (1) 


where R(w) is the frequency response of the media to 

a pulse P(w) and H(w) is the transfer function. It 
has generally been assumed that crust of the Earth has 
been made of a number of layers separated by first and 
second order discontinuities. For such a layered 
medium, transfer function calculations can be performed 
by a number of methods. The Haskell-Thomson (1962) 
Matrix formulation gives a rapid and exact’ calculation 
of a transfer function for a plane parallel multilayered 
mediumin which seismic anisotropy may exist. Gilbert 
and Backus (1966) formulated Sarasa matrix method in 


which the Haskell-Thomson method is a special case. 
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One advantage of this method is that a medium with second 
order discontinuities is allowed (Richards, 1971). 

Ray theory can be used to calculate the response of 
layered media (Hron, Kanasewich, Alpaslan, 1973). 

A ray solution can be applied successfully to dipping 
or curvilinear interfaces provided that a sufficiently 
large number of rays as included in the response. 

In this thesis both the Haskell-Thomson matrix formu- 
lation and a ray theory approach have been used to 
determine the response of multilayered media to a 
plane wave impinging on the bottom. 

In equation (1), the nature of the source pulse 
is generally not known in practice and an incorrect 
assumption can seriously alter the response function. 
If the vertical and radial components of particle 
motion can be calculated, then the ratio of these two 
components will be independent of the source pulse: 


r i 
Equation (2) is called the "transfer ratio" or the 
"spectral ratio". Both terms are used in this .thesis. 
In practice the vertical and horizontal components of 
ground velocity are generally known, then the left 
hand side of equation (2) can be determined and compared 
to the right hand side of the same equation which repre- 


sents a particular model. 
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2. Haskell-Thomson matrix formulation 


Periodic solutions of elastic equations of 


motion for the ee layer may be obtained by combining 


the dilatational wave solution in the x-z plane: 


ae ($3) ey (S) = expliu(t-%)][e exp(-iur, 2/c) 


+ cee sd A alr Ch) 


with the rotational wave solution: 


Ld 3 ; 
OQ =glss- ayl= expliw(t - S)112.4 ei NA 


+ Q ou St eka (2) 


where u,w are the horizontal and vertical radial 


components of displacement, ry: Ye, are che scot’. Of 


B 


the angles of direction of propagation for P or S 


9 2 7) are amplitudes 


aves from the z axis 
bs ag eo oa’) ou 
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of up and down going plane waves and c is the apparent 
(horizontal) phase velocity. 

For . layered media, the boundary 
conditions at the interface between two layers are 
that the displacements u, w, and the normal and shear 


stresses, Paz! Poy are continuous. Haskell (1953) 


has found that the velocity displacement and stress 
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Orientation of axes and numbering of layers 


and interfaces. 
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where the factor exp [iw(t - =) ] is: omitted. U, wih are 
time derivatives of displacements and Da isa 4x4 
matrix whose elements are given explicitly. Equation 


(3) may be rewritten as 
Wt ra Ba | (4) 


where Via and An are the column matrices. At the top of 
the layer, due to the boundary conditions, the equation 


to be satisfied is 
V = Ent (5) 


where matrix En is obtained by substituting z=0 in 
matrix Dat Substitution of equation (5) in equation 


(4) for Hn yields: 


Re = z 
Vin = Dt Vad re eae ol r (6) 


This recursive relation can be repeated from the oe 


layer until the surface is reached: 


Van = An@n-i'** 41% : (7) 


Using equation (5), amplitudes can be expressed as 
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where matrices (a) are known for each layer. If an 
incident plane P wave with unit amplitude is desired 


in the half space, this is obtained by setting em = 1 
n 
and nod. ss oa, = a8 = 0 is the column matrix Hie 


Stresses at the free surface are made to vanish by 

setting the components of V_:(P__) = (P_.) = 0. 
fe) 22% ZX 5 

Then the surface displacement can be calculated in 


terms of the input amplitude for vertical and horizontal 


components: 


2 
- 2c (Jy5 - J35) 
TV, (w) BS 5 
CD 
n 
(9) 
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i An 31 
TW, (w) <= 2 
a, Dry 
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n en-12n-2° °° °F 
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(44 - F51) (339 - Fga) — (Jy 2 ~4Yn9) 31 - Fqy) 


= csci = csci 
Cc On s 1) Ba B 


m m 


(i). is the angle of incidence in nee layer. 


Similarly for an incident plane S wave: 
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TV, (w) = —+2** 
(10) 

2035. = ds) 

Zi BI 


where Ere 2B°/c*, a and g are the P and S wave velo- 
cities respectively. 

When P and S wave velocities, density and the 
thickness of the layers are known, equations (9) and 
(10) can be evaluated numerically to obtain the 
response of a layered media to an input P or S waves. 
The behavior of the transfer function has been inves- 
tigated previously (Alpaslan 1968, Sprenke 1972). 

The modulus of the transfer function has several peaks 
and troughs in the frequency domain whose amplitude and 
peak position change with layer thickness and the phy- 
Sical properties of the medium. The density of the 
layers and the angle of incidence do not cause a 
significant shift of the peaks but the amplitude 
changes can be notable. Therefore, for a crustal 
investigation using the transfer function method, a 
number of seismic sources located at various distances 
Can be utilized. For a given velocity, the density of 
the layer can be approximated closely by using empiri- 
cal data. In this thesis a velocity-density curve 


(Nafe-Drake, 1959) was used to specify the densities 
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of various layers. By varying the parameters of the 
media a good match may be found between experimental 
data and several non-unique theoretical models. 
Within a given frequency band a number of models can 
be found to fit experimental data, This method is 
not a direct method for crustal investigation, however 
it can be applied as a diagnostic tool to choose 
between a number of models. 

A Fourier transform of the transfer function is 
called an "impulse response". Since all phases can 
be considered as a delayed and modified form of the 
first arrival, all phases form a time series in which 


the delay times, tT, introduce an’ oscillation period: 


: (et) 


Therefore large delay times correspond to long period 
waves while short delay times contribute to the high 
frequencies. Path differences between one or more 
phases will determine an interference pattern in the 
frequency domain. For example, if At is the path 
difference between two rays arriving at the same 
receiver, then the phase at a particular frequency (f) 


will be: 


9 = 2n(fAt + n) F 1? = POR 2 eters (12) 


When a number of rays are considered, phase contributions 


are much more complicated. The Haskell-Thomson matrix 
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formulation yields this interference pattern ina 


straightforward manner. 


34 Computation of experimental spectra 


In this thesis, earthquake data recorded by the 
Variable Aperture Seismic Array (VASA) in central 
Alberta, were used to compute the experimental spectra 
and to determine the transmission characteristics of 
the crust in this region. Earthguake data have three 
components of ground motion. One is the vertical and 
others are the horizontal components oriented in the 
N-S and E-W directions. For this study, distant earth- 
quakes were chosen (A>50°), It can be shown that the surface 
of a wave front increases with the square of increasing 
distance from a source. Therefore plane wave assumptions 
are valid for these earthquakes provided that the wave 
front is not diffracted significantly by deeper struc- 
tures of the Earth. Discontinuities in the mantle are 
well separated from each other. The P wave velocity in 
the mantle drops slowly from about 8.0 km/sec, forming 
a low velocity region at a depth of about 100-120 km. 
Thereafter the velocity increases smoothly except at 
depths of 400 and 700 km where there are second order 
discontinuities. According to wave theory, there is 


insignificant conversion of elastic waves at these 
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discontinuities for the frequencies of interest. 
Therefore we can assume that a pure plane P wave is 
incident at the bottom of the crust. Any deviation 
of the wave front from a plane wave will introduce 
errors into the spectral ratio method due to the 
spatial variation introduced. Earthquake sources 
at about 90°- 98° away were also used. According 
to the geometrical ray theory, these events yield 
waves which graze the mantle-core boundary. It 

was assumed that diffraction of elastic waves due 
to diffraction of the mantle-core interface is 
negligible since waves with a wavelength of 4-80 km 
were being used. 

Sixty seconds of digitally recorded data for 
each of three components have been taken for spectral 
analysis. The signal level in every channel has been 
found to vary rapidly. These variations have been 
corrected by bringing signal level down to the mean 
value. The data usually starts one second earlier 
than the first Byrarrival time. Using the thorpzontal 
components (ry (t), rr, (t)) and station-source azimuth 
())*- YTadiat (xr, (t)) and transverse (rp (t) ) components 


were calculated: 


rp (t) ry (t) cos $ + rp (t) sin 9 
(13) 


rp (t) = r,, (t) sin 9 - r,(t) GOS +0. « 
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The beginning and end of the data were tapered by a 
cosine bell in order to reduce an amplitude contribu- 
tion of a sharp discontinuity at all frequencies. 
Using a fast Fourier transform algorithm (Gentlemen 
and Sande, 1966), raw power spectrum was calculated 
from the Fourier coefficients. The finite length of 


data available imposes a time window (W(t)): 
r, (t) = r(t) W(t) (14) 


where the window length can be varied depending on 
various purposes. After 60 seconds of a tapered 
rectangular window, another window was applied to 
the auto-correlation function in the time lag domain 
(t). The window chosen is the Parzen function which 


has a form of 


2 eS 

iggy +6 (ith 3 aie” mn 

Tm Tm 

[| : 1 < 

W(t) = Cie <a) ae lattes Th 
m 
0 isa) es am 
(25) 


where Th. is the maximum time lag. Such a window 
increases the reliability of the spectra and the 
Parzen window has relatively small side lobes and 
yields a positive definite spectrum. Then the real 
Fourier transform of the auto-correlation function 


will give a smoothed spectrum: 
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In this procedure the window can be varied in length to 
eliminate any additional phases such as near source 

and core reflections. But the window length must be 
sufficiently large to obtain the proper resolution. 

For example, the window length used here provides a 
frequency resolution which will separate peaks that 

are 0.025 - 0.033 HZ Apart. ,in vnany events, additional 
phases were included for the sake of resolution. It 

is shown here that by including additional phases, some 
valuable information can be extracted using spectral 
analysis. A Parzen window yields a smooth reliable 
auto-correlation, however it gives poor resolution 

at low frequencies. Therefore, for long-period studies, 
another type of window may be recommended. . 

In equation (1), the transfer function can be 
written in terms of a number of components such as the 
instrumental transfer function, directional transfer 
function etc. (In practice, a knowledge of the instru-~ 
mental transfer function will be of great help in the 
interpretation of the spectrum. However, these func- 
tions have negligible effect on spectral ratios. But 
the gain of the channels must be known for spectral 
ratio studies. Gains can easily be calculated using 
calibration pulses and spectral component can be scaled 


before rotation: 
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(7) 
For the data used in this thesis, there was incomplete 
information available about instrumentation. Therefore 


the amplitudes of the spectral ratios are imprecise. 


&. Efficiency of the spectral ratio method 


As seen in equation (16), the experimental 
spectrum is always obtained by a window and taking 
the ratio of components will not lead to a cancellation 


of this window since a convolution is involved: 


Ran) (w) *W(w) 
Sven ave navn : (18) 


Rep (w) Rp (w) «Wi (w) 


According to this method the response is considered as 


a convolution of one or more functions (equation (1)): 


Roy (w) , [H,(w) +P (w) 1*W(w) 


oO: Sea a . (19> 
Rop t¥ [H, (w) +P (w) ] *W(w) 

As an approximation, if we can write 
[H(w) -P(w) ]*W(w) = [H(w) *W(w) ] -P(w) (20) 


then in equation (19), the unknown source function will 
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cancel, leaving only the crustal transfer ratio function 


where the window effect is not negligible due to 


convolution. Generally, the window in the frequency 
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Figure (i.2 


Theoretical spectral ratios using Parzen window 
and various pulse functions. Pulses effect is 
negligible in the frequency band but minor 
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differences exist towards higher frequencies. 


19a 


8 
8 8 
g PULSE SPECTRUM 
8 
a 
5 SPECTRAL RATIO M258 
8 
3 
8 0. 4.0 6.00 
$ 
¥ 
r8 
x 
: * 
8 60 0.60 1.80 FREQUENCY 3.20 4.00 4.00 s.60 
z 
8 
\ 
a 
s ge eee 
.00 0.20 0.40 0.80 0.60 1.00 1.20 outs) “gio ~C a aT eT CT aC 
8 
a 
8 8 
g SPECTRAL RAIA a, PULSE SPEC1RUM 
2 PULSE 
; i 
‘ =8 
3 24 M=6.3 
Be =8| 
al - 
& log 
8 8 
z cl in a0 
SECONDS “< 
8 8 me 
¢ S00 : oy eer = 
FREQUEAL F 
8 
00 0.20 0.40 0.60 0,80 1.00 1.20 ‘Frequency 100 2.00 2.20 2.40 2.60 2.80 300 3.20 
8 
F 
8 
H SPECTRAL RATIO $, PULSE SPECTRUM 
8 PULSE 
i 
78 
=| 
fa] 
8 88 
ES 53+ 
E E 
8 
| 


8 
.00 0.60 2.40 4.00 4.60 5.60 
FREQUENCY 
1.20 1.40 ro 00 2100 220 2.40 2.60 200. 3.00 9.20 
¥ FREQUENCY 
= ~ 
8 PULSE SPECTRUM 
: PULSE 
“8 
ee 
3 
Be 
a 
Ms 0.00 2.00 6.0 
, SECONDS 
eh, 
4.00 4.00 s$.00 


eeel 


_ Mae 


-~ 


aaa ees . | 
| 
M 


f \ | H 
sok i ii i) | 
—- nl i : ik Ta vie ee ay i 
i, ¥ 4) i 


_~ ‘ iM 
a. tha . <aae a) dca 
j i five +9 ) 
oh TSH | 
> 
My 
Si ual od 
ae BPN ian le Fae oe cep 


mT aa lUF 


Bae 6 
ee ee LS 


ear ee i age $e 


”\ gparoate’ wea 


i 
. {A 
7 eee 
' i] 
, ore s ; ; et | 
ds iA oe | he Vr 7 
So My) Wires ean eas OVO) ECL ht ee 
te Ape ana eat (8) ay OO nT 


if st 


Me Na on Ne be = a 
Pcie an 


eee 


snub ut it A oO alia il 0) ee 
— i ee al Mee rath i 


' ; ‘ 
a) _ Ba 
hs 
»% i Pi ae \ i ela bal Ay Vaan a [ : 
fT 4 ? , y | . 


VR ey aN sae 


en 1 eee | 
oe 7 My “he 


ve a i 


Figure 1.3 


Effect of additional phases on the spectral ratios. 
(In figure b, broken line represents a spectral 


ratio with a phase velocity of 23,6 km/sec.) 
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domain is not a simple function and its effect on the 
ratios cannot be estimated theoretically. Therefore 
to obtain a better comparison between theory and 
experiment, the theoretical models have to be smoothed 
also by a Similar window. 

The approximation in equation (20) can be shown 
to be true for a rectangular time window and the smooth- 
ing of a theoretical model can be achieved by a double 
Fourier technique (Leblanc, 1967). In an earlier work 
(Alpaslan, 1968), this method was applied successfully 
to obtain a theoretical transfer function. However, 
when more complicated windows are used one must be 
assured of convergence and also that no shift in fre- 
quency is introduced. As a first step equation (19) 
was used for a model with a Parzen window and various 
pulse functions as obtained from real earthquake data 
using a homomorphic decaiace lesion technigue. The 
results are displayed in Figure 2. Pulse functions 
used vary significantly in amplitude and frequency 
content. As can be seen spectral ratios for every case 
are similar. The amplitudes and shapes of the ratios 
at frequencies less than 1.40 Hz are exactly the same 
for every case. However, towards the higher frequencies 
small differences exist due to the variable spectral 
content of the source pulses. In general, it can be 


said that variation in the pulse produces negligible 
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variations in the spectral ratio. The most reliable 
ratios were obtained only at low frequencies. [In 
this process it is essential that the window length 
is larger than the pulse duration in time. 

In practice a seismologist can rarely obtain 
data which satisfies all the assumptions discussed 
earlier. For example, after the P arrival, additional 
phases such as pP and PcP arrive within the window 
length. These multiple pulses create holes and dis- 
tortions in the experimental spectra destroying valua- 
ble information. If one can choose earthquakes hav- 
ing a great depth (>150 km) and distance (A>55°), 
about 40-50 seconds of data will be free of multiple 
phases. In this thesis if such data were chosen, only 
one or two events would be available. Therefore these 
phases were included in this study, but many spectral 
ratios were combined to obtain ratio that approached 
the true one. In Figure 2, many models with additional 
phases have been displayed. Additional phases were 
introduced to the theoretical Haskell-Thomson formula- 
1 Ons 
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R(w) = Re (w) x Re, Wwe Pre, (ule eo ico 


At is the time delay of a particular phase and Cy 1Co 13 


are the phase velocities. All pulses used here were 
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taken from real data using a tics SUG er ae technique, 
therefore amplitude of these phases were not assumed. 
In Figure 3.b, all phases are assumed to have the same 
phase velocity as it can be seen that the spectral 
ratios are only affected in amplitude by comparison 
with a single pulse case (Figure 3.a). In Figure 3.c,d,e, 
the phase velocity of pP is increased steadily. It is 
clear that the spectral ratios change their appearance 
due to the increasing difference between main P and 
pP phase velocities and the delay time will determine 
amount of distortion. Therefore, it is concluded that 
the spectral ratios are still useful even when multiple 
pulses exist within the window length provided that the 
phase velocity of pulses does not differ significantly. 
For spectral ratio analysis, data with high 
eat to noise ratio must be chosen. ree data will 
produce extra peaks in the spectral ratio which cannot 
be modelled by a structural effect. This is one of 
the possible reasons one obtains different spectral 


ratios for the same station when using various events. 
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Ite 4 ANLSOtropy 
Toe Lnerouuce Lon 


A large amount of evidence about the seismic 
anisotropy of certain regions of the solid Earth have 
been accumulated for a long time. Neuman (1930), 
during a study of various earthquakes, observed an 
impulsive activity before the main S arrival. He com- 
pared this motion to the Love wave type and after the 
main S arrival, he identified Rayleigh wave motion. 
Byerly (1934-38) has noted the change of S motion from 
transverse to the vertical plane. He suggested that a 
double refraction phenomenon occurs in the surface 
layers of the Earth, Measurements of seismic anisotropy 
on rock samples confirmed the variation of the velocity 
with direction (Koefford et al, 1963). The degree of 
anisotropy was found to vary for a different rock sam- 
ples. Khalevin and Koshkina (1966) measured the sound 
velocity in rock samples collected from the Ural region 
of the U.S.S.R. Among Pease Voor samples, metamorphic 
rons showed 300 % anisotropy. The variation of velocity 
is observed in three directions; one is the perpendicular 
to the cleavage; the second is parallel to the cleavage 
and parallel to the lineation; and the third direction 
is parallel to the cleavage and perpendicular to the 


lineation. Their results showed that the velocity in 
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the second direction is larger than that in other 
directions. Compact textured samples (granoblastic) 
showed small anisotropy. 

Small scale anisotropy has been measured by a 
number of authors in situ (McCallum (1932), Wheatherby 
(1934), Ricker (1953), Vanderstoep (1966), Cherry et 
al (1968)). Yegorkino (1969) used recordings of a 
distant earthquake at four stations separated by 7-10 
km, to determine the anisotropy in a sedimentary 
layer in the northern region of Caspian Basin of the 
U.S.S.R., using travel time and polarization analysis 
for converted phases (PS..-PS_) at the bottom of the 
sedimentary layer. He was able to measure 0.5-0.7 
seconds time differences in arrival times of PS, and 


V 


y Wave velocity was 


wave velocity in this 


H 
10-15 % larger~ than the S 


PS., waves. He concluded that the § 
V 
sedimentary layer. 
Anisotropy has been shown as a possible cause 
for obtaining "non geological" boundaries in a number 
of seismic field measurements (Layat et al, 1961). 
Khalevin and Tavrin (1965) could not find any corre- 
lation between boundaries obtained with the deep 
seismic sounding method and those established by 
stratigraphic and lithographic studies in the upper 
part of the crust (450-850 m) in Urals of the U.S.S.R. 


They believed that the main reasons for this was the 
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varying degree of metamorphism in the rocks of the 
upper part of the crust. 

A theoretical description of general anisotropy 
can be made in terms of 21 elastic parameters (Love, 
1944). Backus (1970) offered an interpretation of 
these parameters. However, a complete physical inter- 
pretation of all 21 parameters is not possible at this 
time. Using various symmetry relations, a number of 
Simpler forms of anisotropy can be stated. The simplest 
form of anisotropy allows a physical interpretation and 
might partially explain anisotropy measurements. One 
such type of anisotropy is the transverse isotropy in 
which the propagation of elastic waves perpendicular 
to an axis of symmetry is equivalent to waves propagat- 
ing in an isotropic media. This symmetry axis is called 
the unique axis. Only five elastic parameters are 
sufficient to describe this kind of anisotropy in which 
SH and SV waves can travel with different velocities 
and P waves can have three different velocities of 
propagation (Stonely, 1949). Surface waves exist in 
a transversely isotropic media as is the isotropic 
case (Buchwald, 1961). A number of features of surface 
wave dispersion have been explained by transverse 
isotropy (Kovach, 1965; Kanimura, 1966; Crampin, 1966). 

Transverse isotropy can generally be divided 


into two groups: (a) transverse isotropy with a vertical 
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unique axis. In this group SH and SV waves travel with 
a different velocity. (b) transverse isotropy with a 
symmetry axis other than vertical. This kind of 
anisotropy predicts azimuthal variation in P wave 
velocity. Hess (1964) reported Po velocity variation 
with azimuth in the Mendocino fractured zone, off 
California. Raitt (1969) measured anisotropy as a 
function of azimuth in the Pacific ocean. Backus (1965) 
derived a relation for an anisotropic medium in which the 
P wave velocity varies with azimuth. He showed that 
measurements of velocity in five different azimuths 
will determine the anisotropy. This method has been 
applied by Christensen, Crosson (1968) and Raitt (1969). 
It can be stated that the formation of transverse 
isotropy in the real Barth may form in two ways. One 
is the sedimentation and metamorphism of Earth material 
(Khalevin et al, 1965). The other is the orientation 
of crystals under anisotropic stresses. Macdonald (1960) 
studied orientation of anisotropic minerals in a stress 
field. He showed that in a stress field each mineral 
tries to take an orientation in which the potential 
energy of the system is minimum. The direction of 
orientation is shown to be the same as the one with 
the lower elastic wave velocity in the minerals. 


Christensen (1965) measured the compressional wave 
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velocity in various metamorphic rocks at high pressures. 
He observed a variation of velocity with direction as 

a result of mineral orientation. Silaeva and Bayuk 
(1967) measured the P wave velocity and anisotropy in 
rock samples at a high pressure and temperatures. They 
have concluded that anisotropy depends on the stress 
state of the rocks. Corlet and Emery (1959) measured 
the vertical and horizontal pressure in sedimentary 
rocks. They’ have found that the horizontal pressure 

is 0.5-0.75 of the vertical pressure . In tectonically 
young and active regions comparatively low sub-crustal 
velocities of longitudinal waves (7.5-7.9 km/s) are 
partly due to the incomplete orientation of the medium 
(Babuska, 1970). Generally, in tectonically active 
regions such as young mountain ranges, suboceanic 
trenches etc. at least triaxial pressure a proposed. 
Therefore one may have a transverse isotropy with a 
unique axis not vertical. For a tectonically quiet 
structure, assuming a vertical biaxial pressure, 
transverse isotropy with a vertical axis is the usual 
form of anisotropy. 

In this thesis the possible implications of 
anisotropy on the transfer function have been inves- 
tigated. The purpose of this study is not to determine 
the anisotropy in an Alberta crustal section but to 


clarify the uncertainty imposed by the assumption of 
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isotropy during the determination of crustal parameters 
using a transfer ratio method. The study of anisotropy 
has been achieved in two steps. In the first step, 
experimental transfer functions calculated at each 
station in Alberta have been classified according to 

the direction of wave approach. Spectral analysis 

has been carried out for each quadrant at every station. 
In this way the possible structural anisotropy in 
azimuth is investigated. The result of this study 
depends on the number of reliable experimental data. 

The conclusions and diagrams for this study are pre- 
sented in another section where interpretation of 
experimental results are discussed. In the second 

step which is presented in this section, the transfer 
ratios for transversely isotropic media with a vertical- 
unique axis have been calculated using a matrix method. 
Models are assumed to be two-dimensional, homogeneous, 
plane and parallel layered media. Various degrees of 
anisotropic models have been compared with isotropic 
models and differences have been discussed in terms of 


the structural interpretation. 


Ze Computation of a plane wave solution for a trans- 
versely isotropic multilayered media 
The equation of motion for elastic waves may be 


written (neglecting body forces) as 
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xu, oP. . 
fe) ) = 5%. J = 1,273 €is) 
ot J 


where p is density, U is displacement. We assume an 
elastic body that is, it produces internal stresses when 
subjected to an external force that causes deformation. 


The stress components Pi can be written as a function of 


strain (e._) 


Pq 
Pia = PS 2 
is i3 ‘pq? (2) 
where 
1 dU 0U 
og 2 kes + a eee es 


If the strain components are small and we can assume 
absolute elasticity, which means that the stress 
becomes zero in the absence of deformation, then 
expansion of equation (2) into a Maclaurin series 
leads to a linear relationship between stress and 


strain, which is known as Hooke's law 


Pay = Pi ipaene Ly Dre pee (3) 
where 
Oils 
rns = Fses are called elastic parameters. 
1jpq one) 
PT e +0 


se 


P..ande are second rank tensors and c.. LS a 
a pq 13~9 


three-dimensional fourth-rank tensor, independent of 


the choice of a cartesian coordinate system. This has 


.srstemezag Sitesis belts 92s 


oe 


no omvees SW ,. 


.soitsmoteat aseimo teas | sna 3 


m > 


Pai) 


to anigonn? 5 Be esti ad. aaa L 


= net 


(S) 


euetts silt jas een siy tee “4 


ine 


itedy , noisenaoteb ‘to bers 


bas eeexte neewded otal ey 
é wel at esook ae ‘i 


‘ 


¢f) ES, h = aA: 


81 independent components. Then the equation of motion 


becomes 
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The 81 parameters are reduced to 21 by the symmetry 


relations 
Cues = C.; =O. Se Clon: : (5) 
1)P9 J1Pq 1390p pgqi)j 
The elastic parameters Spiess are dependent on whether 


the process is adiabatic or isothermal. It is 
necessary to assume a positive definite function ¢ 


exists (Love 1944), such that 
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where 


d= i Cc e..e 
Ziaesgpa 49.4pq 


The function © is termed as elastic potential per unit 
volume. It is the density of potential energy in an 
elastically deformed body and is defined in terms of 
the thermodynamical conditions of the body. 

A plane wave solution of equation (4) can be 


written as 


U, = a,expliw(t - MX)! (7) 
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: = (Ak, , W is the radial frequency, k is 
the propagation vector of plane waves in a cartesian 
coordinate system. Then the equation of motion (4) 


will be 


[= eee o,f eer io of PON sy 
Equating the determinant of equation (8) to zero, the 
displacement components can be calculated. For 21 
elastic parameters, the evaluation of these equations 
is tedious. However, for a particular crystal system, 
the determinant can easily be expanded. 

Transverse isotropy has a symmetry which is 
equivalent to a hexagonal crystal. A six fold symmetry 
axis exists. The elastic parameters are unchanged by 
rotation about this axis and nee of the parameters 
vanish or become equal to each other, leaving only five 
independent elastic modulii. In a cartesian coordinate 
system where the x, axis is downward, the components 


of the elastic parameters will be 
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A separation of the equation of motion into 
longitudinal and transverse components is not possi- 
ble. Therefore, the distinction between P and S waves, 
cannot, in general, be done, except in certain direc- 
tions. An expansion of the determinant of equation 
(8) yields relations between the velocity of wave 
propagation and the elastic moduli (Stoneley, 1949). 
Considering wave propagation in certain directions, an 
interpretation of five elastic parameters can be made. 
Following Stoneley (1949) and Anderson (1961), some 
directions of wave propagation in a transversely iso- 
tropic media are comparable to the isotropic case: 

a) For a direction of propagation along the unique 
axis (x3), pure compressional and shear waves with 
horizontal particle motion are the solutions of the 


determinant of the coefficients: 


al al 
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boy | 2325 
P = SV = ene, 
Uera (res ane. Gevarm (10) 
b) Propagation in any direction perpendicular to 


the unique axis, yields the following solutions: 
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PV, PH can be considered as P wave velocities in the 
vertical and horizontal directions. SV is the 
velocity of a vertically travelling S wave with hori- 


zontal particle motion and SH,, corresponds to an S 


H 


wave travelling perpendicular to the direction of SV. 
waves with a horizontal particle motion. The fifth 
parameter can be-determined if the velocity in any 
other direction is known (usually 45° from the vertical 
axis): 


apo eye ie (12) 
1133 Orn 2328 


where Py is the P wave velocity in a direction other 


than the ones parallel to the axis, and 
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For a real velocity, the determinant of equation 
(8) yields six roots for every These roots might 
be real or form complex conjugate pairs. Every root 
contributes to the displacement (equation (7)) and the 
theoretical implications ,;of these ywroots have been 
investigated by Synge (1957), Kraut (1963) and by a 
number of other authors. However, the solution of a 


wave equation can, in general, be written in terms of 
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existing up and down travelling waves (Anderson, 1961) 


in the (X)-X3) plane: 


kr X3 ~krx3 hrax -kr_x 
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(23) 


where the common factor exp[i(wt - kx, )] is omitted. 
Each term in equation (13) represents plane waves 


travelling in x, direction which is pointing downwards. 


=) 


U U are the amplitudes of up and down 
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going P and S waves. lp and r. are the cot. of angle 


up’ 


of direction of) P wand®’S waves, (from «x. ‘direction. + k = w/c, 
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and c is the phase velocity. 

By substitution of equation (13) into equation 
of motion and equation (6), displacement and stress 
components can be written in terms of input amplitudes. 
It is more convenient to use dimensionless velocity for 
the displacements (Harkrider, Anderson, 1964). Then at 


the meh interface the solution can be written 
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+1 
= U + (-1) U 
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and 
(c + c RAS 
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The stress components will be 


em ik [ix,U.jcos P-X,U,5sin P+ iX,U,3cos Q- 


- XU. ,sin ol 


Gry 
(Pi 3), = ik[iY,U.,sinP- YpUogcos P + iY,U,3sin Qt 
+ YoU. ,cos Gh 
where 
A Shes 3c lta ee eerie 
h = P or S 
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For a plane parallel multi-layered media, equations 
(14) and (15) can be written in matrix form. Recur- 
Sive relation for a number of interfaces can be 

formed in the same way as it has been done in section 
L, to £ind the surface amplitude in terms of the input 
amplitude and the vertical and horizontal transfer 
functions can similarly be written for a (n-l) layered 


th 


media, the n layer being a half-space (m = 0,1,2,..n). 


Then, the transfer function for an incident P wave 
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D and matrix J are similar to the ones calculated in 


section I. 


3, (transfer Tatios for transversely isotropic layered 
media 

Equations (16) and (17) are the surface responses 
of a transversely isotropic multi-layered media to an 
input P wave with unit amplitude in a form convenient 
for numerical evaluation for a number of layers. In 
Figure 1 the transfer ratios which are the ratios of 
the vertical component to the horizontal one, have been 
shown for an Alberta crustal model using two different 
angles of incidence. Anisotropic factors have been 


shown as the ratio of velocities: 
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Figure Ee 


Effect of transverse anisotropy on transfer 
ratio as a function of angle of incidence: 
PHT =.-0'..83), (BTA = 0.93 an thes crust. 


Broken line indicates isotropic model. 
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Figure, £i.2 


Effect of transverse anisotropy on transfer 
ratio as a function of angle of incidence: 
PHI .= 0.83, ETA = 0.91 in the crust. 


Broken line indicates isotropic model. 
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All models presented in this thesis are two dimensional, 
therefore the factor (XI) should not be considered on 
these figures. Models presented in Figure 1 have 10 % 


anisotropy in which the P,, velocity is higher than P 


H V 


velocity and the half space is assumed to be isotropic. 
By comparison with the isotropic case, peaks and troughs 
are shifted to higher frequencies indicating shallower 
or higher velocity structure. This can be explained 

in terms of energy propagation, since in an anisotropic 
media, energy does not necessarily propagate perpendi- 
cular to the wavefront. It may even propagate parallel 
to the wavefront with a higher velocity (Lighthill, 1960). 
If the energy flux vector or Poynting vector is divided 
by the total energy density of the wave, the velocity 

of energy apepaGoton (group velocity) is obtained. 

The relation between group velocity and phase velocity 


is found to be (Musgrave, 1954): 


the n = Vp | (18) 


where Vg is* group “velocity, Vp is the phase velocity 
and n is the normal to the wavefront. In an aniso- 
tropic media the direction of displacement vectors 
and the Poynting vector must be considered at the 
same time when reflection and transmission of an 


elastic wave is studied. The critical angle phenomenon 
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is more complex in an anisotropic media and needs 
more explanation in terms of these vectors. 

It is clear that when the horizontal distance 
travelled by the wave in an anisotropic media, is 
Significant with respect to a vertically travell- 
ing wave, the effect of anisotropy will be greater, 
when transmission of a plahe wave through an multi- 
layered anisotropic media is considered. The angle 
of wave approach at the bottom of the media is 
important since a variation of velocity with direction 
will be felt easily by the wave for grazing angles of 
incidence. A larger number of models, than shown in 
this section, have been calculated. Even small aniso- 
tropy introduces important changes in the shape of the 
transfer ratios for large angles of incidence. How- 
ever, for experimental transfer ratio studies, distant 
earthquakes are chosen in this thesis. Most of these 
teleseismic events should yield waves with an angle of 
incidence less than 35 degrees at the bottom of Alberta 
crust, according to the Jeftreysnul ten tables (Phinney, 
1964). 

The anisotropy factor (n) will determine the 
state of transverse isotropy in a model. Models 
shown in Figure 1 have an intermediate velocity (P,) 
which is equal to the vertical velocity (n = 1). 


This sharp velocity change in horizontal direction 
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Figure if23 - 


Influence of upper mantle anisotropy on 
the wave transmission. Broken line indicates 


isotropic model. 
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corresponds to highly laminated media. The degree of 
lamination can be introduced by a factor n. For 
example, models shown in Figure 2 have a smooth velocity 
variation from the horizontal to the vertical direction 
(7 = 0.91). This model showed no significant shift of 
peaks in frequency but only amplitude changes. From 

the comparison of Figure 1 with Figure 2, it can be 
stated that the teaneter ratios are much more influenced 
when n = 1. The transverse isotropy introduces no addi- 
tional peaks into the transfer ratios by comparison with 
the isotropic case. The effect of transverse isotropy 
on transfer ratios is not so severe for low frequencies. 
For frequencies less than 0.26 Hz, the maximum frequency 
shift in peak position is observed to be 3 mHz and the 
amplitude changes are found to be 20-30 % from the 
isotropic case. 

Since the angle of wave approach at the bottom of 
themedia is important for transversely isotropic media, 
the anisotropy of the half space will be one of the 
contributing factors. In Figure 3, an’ isotropic Crus— 
tal model over the transversely isotropic half space 
has been studied. This transfer ratio has been calcu- 
lated assuming that a purely incident P wave exists. 
Existence of a purely P wave in an anisotropic media 
may be questionable, however, it can be assumed that 


due to an anisotropic process at the source, deep 
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Figure I1l.4 


Comparison of an experimental spectral ratio 
with an Alberta crustal model (ALTACRT-1). 


For unsmoothed models amplitudes refer to 


Square root of power. 
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isotropic structure of the Earth or by some kind of 
mechanism, a pure P wave is introduced into a trans- 
versely isotropic space below the crust. Investiga- 
tion of these assumptions will not be done here since 
this model is only of theoretical interest at this 
time. This showed that 50 % amplitude changes (Figure 
3) with respect to the isotropic wave indicating the 
importance of the angle of incidence. 

For long period elastic waves, if the amplitude 
changes can be ignored, the effect of transverse iso- 
tropy on the transfer ratios is negligible. However, 
for short period waves, a study of transfer ratios 
is no longer independent of the angle of incidence 
and anisotropy may be one of the significant factors 
to overcome some difficulties in application of this 
method to the real data (Ellis and Basham, 1968). 

In Figure 4.b-c, an Alberta crustal model (Cumming 

and Kanasewich, 1966) has been compared with an experi- 
mental transfer ratio. The experimental spectral 

ratios have been calculated using an earthquake recorded 
at the EDM Station in Alberta (section I). This 
spectral ratio is smoothed by a 30 second Parzen window 
in order to increase the reliability of the spectrum. 

By comparison with the isotropic model, it can be said 
that the experimental and theoretical transfer ratios 


are not completely similar. However, the anisotropic 
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model makes a comparison much easier. This anisotropic 
crustal model has a vertical Poisson's ratio oO = 0.24 


and a horizontal one On = 0.35 in the crustal layers. 
There is no other evidence for various Poisson's ratio 
values in this crustal model, therefore, this model 

must be verified. Both models yield a very poor com- 
parison with the experimental transfer ratio at both 
frequencies. This may mean that the model has to be 
modified. The models shown in Figure 4.b-c were not 
smoothed because small peaks in the experimental curve 
have very low reliability. In Figure 4-a, another model 
with six layers over a half space (Table 1) has been 
displayed. This model is smoothed by a similar window. 
Comparison with the eee nent data is more favorable 
than for a four layer model. This may suggest that 
significant alteration has to be made in the sedimen- 
tary layers of the model. In the Alberta plains, a 
thick sedimentary layer is a common feature of the crust 
and anisotropy in these layers may significantly affect 


the transmission characteristics of the crust. 
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Table 1 


Crustal Model ALTACRT - 2 


h (Thickness) P (km/sec) S (km/sec) ep (gxr/cc) 
LSS 3.00 ISS 2 a0 
0,6 Dez 350 200 
0.5 4.30 2.48 260 
10.0 6cL0 SoZ oe 10 
Zou) 6550 SS Pa /S: 
10:0 Te tS 4.13 +20 
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Pet peeelec’ OL Dipping Layers on Wave Transmission 


ges Introduction 


In the study of the structure of the crust one of 
the simplifying assumptions has always been the existence 
of plane parallel layers. It is most likely that this 
assumption is not satisfied in many crustal sections. 

As far as we know, Kuhn (1961) was the first to study 
the time response of a dipping layered structure using 

a laboratory model. Kane (1966) applied ray methods 

to determine the response of a single dipping layer over 
a half space and Fernandez (1967) discussed the effect 
of dipping layers on the transfer function. Ishii and 
Ellis (1970) formulated the solution of time transformed 
wave equation for a dipping layer over a half space in 
the cylindrical coordinate system. Rogers and Kisslinger 
(1972) have used this ray method to calculate the 
response of a single layered dipping model in the fre- 
quency domain and attempted to determine a dipping 
crustal structure in U.S. They have also made a model 
study. Tsukuda (1972) has attempted to determine the 
response characteristics of a single dipping layered 
analog model. He has concluded that for dip angles 

less than 5-10 degrees, there is no significant change 
in the transfer function for plane P waves. His results 
are concordant with the ones obtained by Rogers and 


Kisslinger. 
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An analog model study can be considered as a 
complete wave solution for that particular media. 
However, an analytic complete wave solution for a 
dipping layered media cannot be formulated as in the 
case of a parallel layered media (section I). In this 
section, a ray-theory approach has been used to evalu- 
ate the response of a multi-dipping layered media in 
the time domain. The transfer ratios have been calcu- 
lated from the Fourier transform of the impulse response. 
The have solution can be expanded into an infinite 
number of series, and each term represents a physical 
ray (Spencer,1960; Cisternas et al, 1973). For numerical 
evaluation, the ray expansion has to be terminated; 
therefore, a comparison between the complete wave and 
a ray solution has to be made to determine the proper 
number of terms in the ray expansion for a suitable 
approximation. The Haskell-Thomson matrix formulation 
(section I) is an exact wave solution and has been used 
to determine the accuracy of a ray expansion. A com- 
parison of the two solutions has been carried out in 
both the time and the frequency domain. After deter- 
mining the optimum number of terms in the ray expansion 
to obtain the best approximate solution economically, 

a ray method has been used to calculate the response 
of a dipping layered media to a plane wave impinging 


on the bottom. Using an efficient algorithm for ray 
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generation (Hron, 1972), the transmission of a plane 
wave through a multi-dipping layered media has also 
been studied. In this section only the transfer func- 
tion studies have been presented as the time domain 


studies will be discussed later. 


ee ey ray-theory approach For computation of the 
response of a layered media 


In an isotropic homogeneous media, the displace- 
ment vector can be written as a sum of longitudinal and 
transverse components (Bullen, 1963). Each component 
satisfies a wave equation of type (section I) 


A inte ee 
AUS = 1) ees 0 (1) 


where A is Laplacian and U, is time derivative of dis- 


t 
placement U(x,y,z,t). A general solution of the wave 
equation for a disturbance travelling in the positive 
direction with a velocity, c, and whose position in 


space is given by a surface W(x,y,z), can be written 


as 


D(x,y,z,t) = H(x,y,z,t) f(t — Veeye2)) (2) 


where the function f(t) represents the disturbance and 
it may be discontinuous at t = 0. The position of the 


surface, W(XiV;z), at any time: (Ct). is given by 


t = (x,y,z) = W(x,yr2) (3) 
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Figure vila. 


Approximation of wave solution by ray expansion. 
The long period spectrum was calculated using a 


time series of 35 seconds long. 
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The expansion of equation (2) in a Taylor series near 
the surface, W(x,y,z), (t = ttAt) leads to the follow- 


ing equation: 


co n 
Ulx.Ve7Z ott) Soepe), a (oe »(At)™ Pit stole, ¥, 2) } 
MeO nay sete 


t=t (x,y,Z) 

(4) 
In this equation the amplitude of the disturbance is 
no longer time dependent. For the particular case of 
a harmonic wave solution, where the amplitude is 
expanded as a power series in (l/iw), it can be taken 


as (Hron and Kanasewich, 1971) 


iw (t-T) 
(At)” £(t-t) = <_—___ ' (5) 
(iw) 


= 
n 


Then the solution can be written in terms of only 
space dependent amplitudes A(x,y,z): 
ok hs eiw (t-T) 
Ux Voz Se) ORC NES Peemper rcter 0 (6) 
n=0 (iw) 
Equation (6) can be generalized for an arbitrary source 
function S,, ft) 
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S(w) is the Fourier transform of S(t) and defined to 
be zero from w= 0 to Woe Equation (7) is known as an 
asymptotic ray expansion. 

If we assume a solution of the wave equation in 
the form of equation (2), for a specific wave whose 
velocity of propagation is v, then substitution of 
equation (2) into equation (1) leads to the following 
equation: 

GHOy any POPS $e (9) 
where it is assumed that velocity is independent of 
space coordinates. Equation (9) is known as the 
"Eikonal" or "characteristic equation" of the wave 
equation and it is a first order differential equation, 


the solutions are 


W(X,Y,Z) = const . C0") 


These are surfaces in three dimensional space, which 

are called wavefronts. Since rays are defined as 

normals to the wavefront, then equation (9) leads to 

the concept of a "ray". Solution of equation (9) must 
also be solution of the wave equation. When it is not, 
some limitations should be imposed on the solution 

such as the wavelength being smaller than layer thick- 
ness, negligible change in the curvature of the wavefront 


and amplitude of the wave over wavelength. In equation 
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(7), terms of (l1/iw) are the corrections to the first 
term due to such limitations. For a plane wave, a 
solution of the Eikonal equation is also a solution 

of the wave equation, provided that wavelengths are 
much smaller than layer thickness. Only n = 0 in 
equation (7) will be sufficient to obtain a wave 
solution using first order ray theory. Proof of this 
is given and discussed by Cerveny and Ravindra (1971). 
Equation (7) can be applied to any media if the source 
function is known at time t = 0. Rays can be drawn in 
Space using Fermat's principle. When there is a dis- 
continuity in parameters of the media, continuity of 
displacement and stress components are the boundary 
condition to be satisfied and the amplitude of the ray 
is predicted by coefficients of reflection and trans- 
mission (Knott, 1899). In an isotropic media it can be 
assumed that energy flow is in the direction of the ray. 
For a plane wave, the cross section of a ray bundle 
(ray tube) is constant at any time and it can easily 

be shown that energy is conserved at any instant. 
Therefore, for this simple case, a geometric expansion 
factor is not necessary (Cerveny and Ravindra, 1971). 
For a ray which is reflected and transmitted a number 
of times, the amplitude, A(x,y,z), in equation (7) will 
be the product of reflection and transmission coeffi- 


cients (Spencer, 1960). For a plane wave solution, 
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the density of rays will be important and the sum of 
existing rays will constitute a solution of the wave 
equation. However, for numerical evaluation of a 
solution, the number of rays to be summed must be 
chosen from an infinite number of rays. 

A program developed by Hron et al (1973) is used to 
determine the optimum number of rays necessary to evaluate 
economically the response of a layered-media to a plane 
wave impinging on the bottom. The validity of the 
solution has been tested by comparison with the one 
calculated using a matrix method (Haskell, 1953) which 
can be considered as a complete wave solution. Rays 
with up to 13 segments for a four layered Alberta 
crustal model (Table 1), have been shown in Table 2. 
Segment number 5 represents the ray impinging on the 
bottom of the layered media. Each segment is repre- 
sented by a number corresponding to the layering 
sequences (Figure 1). For example, ray number 1 has 
five segments and represents the ray transmitted 
directly from bottom to the surface. Rays with a 
large number of segments can be constructed from ray 
number 1, by adding two more segments to it in any 
layer. For this model, the total segment number in 
a ray will be an odd integer starting from 5. The 


ray Classification is based on the number of segments 
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in a ray (NSEG). The maximum number of conversions 

of a ray will be JCONV= NSEG-1. A list of converted 
phases constructed from Table 1 has been shown in 
Table 3. For example, if one conversion is to be 
allowed for ray number 1 which has 5 segments, then 
possible conversions may be at one of the four inter- 
faces: from bottom tothe top’ (JC=1,2,3,4). The total 
number of converted rays with NSEG segments and JCONV 


conversions is given by the equation (Hron, 1971): 


ale mere ieee islet Sd am 
JCONV JCONV! (NSEG -1- JCONV)! ° 


Therefore, the number of rays increases rapidly. For 
example, from the unconverted rays with 13 segments, 
One. can CONStruct about half a million rays. For ia 
plane wave solution, kinematic and dynamic analogs of 
rays are not necessary since phases are well separated 
in time. The number of conversions to be used on each 
set of rays for a satisfactory seismogram is found 


empirically as given in the following table: 


NSEG Seu No. of conversions 
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Generally, ray amplitudes decrease with increasing 

number of reflections as well as number of conversions. 
In this thesis, the minimum amplitude for a ray to be 
included was taken to be 1/300 of the existing maximum 
ray amplitude. The impinging point and offsets on the 
bottom of the layered media has been calculated and the 
only rays of any impinging within 1000 km of the receiver 
have been taken into consideration for a flat layered 
model (Figure 1). For a dipping layered model the 1000 
km region on both sides of the receiver has been reduced 
further to exclude the vertex. The radiation from the wedge 
shaped structure is ignored as it is assumed that waves 
propagate out of the open end of the layer. Using the 
tests given above, accepted rays are determined with 
respect to existing rays. When less than a certain given 
percentage (7 %) of the rays are accepted inisa group, 

no other rays ina higher set are taken into considera- 
tion in order to save computing time. 

The response of a model (Table 1) calculated 
using ray theory has been compared to the one calculated 
by a matrix method in both the time and frequency domain. 
The time domain comparison is presented in another 
chapter and gives a very good comparison over the 
initial interval of the coda, however, small differences 
exist at larger times. The lower end of the frequency 


is the most imprecise because this corresponds to large 
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reverberation times. In Figure 1, the transfer ratio 
for an Alberta crustal model (ALTACRT-1) calculated 
using ray theory has been compared to the one calcu- 
lated using a matrix method. The worst part of this 
comparison is at about 0.04 Hz. In this model, wave- 
lengths are about 6-15 times larger than layer thickness, 
however, there is no change in the position of major 


peak and only amplitude differences exist. 


Table 1 


Alberta Crustal Model ALTACRT-3 
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reverberation times. In Figure 1, the transfer ratio 

for an Alberta crustal model (ALTACRT-1) calculated 
using ray theory has been compared to the one calcu- 
lated using a matrix method. The worst part of this 
comparison is at about 0.04 Hz. In this model, wave- 
lengths are about 6-15 times larger than layer thickness, 
however, there is no change in the position of major 


peak and only amplitude differences exist. 
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3. Dipping layered media and transfer ratios 

A ten degree dip was introduced into model 
ALTACRT-1 into the basal interface (Figure 2). The 
Sign convention was that a positive (negative) sign 
indicated an up-dipping (down-dipping) interface with 
respect to a wave approaching at the bottom of the 
layered media. The difference between a flat anda 
dipping model mainly occurred in the amplitude of the 
transfer ratio and there were no significant frequency 
shifts. In this model a peak at about 1.28 Hz has been 
affected most and the down-dipping model shows a 
decrease in amplitude with respect to the flat layered 
model, which is opposite to amplitude change for an 
up-dipping case. For frequencies larger than 15 Hz, 
very significant differencesappeared. The effect of 
dipping layers on an incident S wave is more serious 
(Figure 3) and this can be explained in terms of 
particle motion. For P waves interfaces dipping 
as much as 10° introduced no significant shift of peaks 
in frequency. However, when low velocity sediments 
have dips about 10°, the most significant effects are 
observed. The important shift of peaks occurs at all 
frequencies as well as amplitude changes. Here, an 
important feature of program must be mentioned; when a 
thin layer at the top of the model is dipping severely, 


the region on which rays impinging will be reduced to 
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a very small value, for example, for +10°, it will be 
(Figure. J) =56.7<x<i1.3 km. Therefore; a number 

of the rays will be excluded in the response. For 
example, in Figure 4, the peaks around 0.90 Hz does 

not exist. When the same layer has a dip of -10°, the 
same important phases such as P.P,P3P5Py do not reach 
to the surface due to the above reasons, raising some 
question about the solution. han all the layers in 

a model are dipping concordant, the interpretation of 
the transfer ratio in terms of real Earth is difficult 
(Figure 4, bottom). For this case a frequency shift 
of about 24 MHz was observed at very low frequencies 
(0.1 Hz) (Figure 5). Generally, about 30-40 seconds 

of experimental data are used to calculate spectral 
ratios, this time length corresponds to a plane wave 
impinging on about 100 km region of the bottom of the 
layered media. Dipping of sedimentary layers as much 
as 10°, brings serious restrictions on the ray solution, 
or this may mean that radiation from the wedge shaped 
structure is not negligible. For this model on which 
sedimentary layer has a severe dip, the peaks at very 
low frequency band were shifted to the higher frequen- 
cies, however, this shift is found to be unstable 

with changing angle of incidence. If Figure 6-d, model 
(i = 35°) showed a shift of peaks towards the low fre- 


quencies with respect to flat layered model, the 
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Figure I1t.2 


Comparison of spectral ratios for a dipping 


model with flat one. 
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interpretation of this amount of shift in terms of 
real Earth can be done if the total thickness of the 


crust is varied about 5 km. 


4. Conclusions 


For dipping interfaces of thick layers with 
a high velocity, the transfer ratios are not affected 
Significantly for dip angles of as much as 10 degrees 
in the frequency range of interest. In flat layered 
media, a layer with a thickness less than one wave- 
length introduces no significant frequency shift into 
the transfer ratio. However, there are significant 
shifts when steeply dipping low velocity layers are 
present. For example, an up-dipping sedimentary layer 
as much as 10 degrees causes a shift to the higher 
frequencies. This is expected, because the low velocity 
layer affects travel time of the wave significantly. 
For a down-dipping sedimentary layer it is expected to 
have a frequency shift to lower frequencies due to 
larger travel time of wave in this layer. However, 
for large down dipping thin layers, the region on which 
wave impinges on the bottom of the layered media is 
very small, therefore the number of rays used is very 
low and the solution cannot be evaluated using this 
program. In both cases, the amount of shift introduced 


by low velocity layer on transfer ratios is not certain 
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Figure III.4 


Effect of steeply dipping sedimentary layer 


on transmission of P waves. 
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Figure III.5 


Effect of up-dipping layers on long period 


transfer ratios for P waves (model ALTACRT-1). 
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Figure III.6 


Effect of dipping interfaces’on long period 
transfer ratios and its variation with angle 


of incidence. 


69a 


8 

3 

es 

§ 

8 

8 

‘03S Ol+ —SE=I Y34SNUUL . 

8 

3 
TSGOW LV 14 

(p) 

se one we wre a ot oe) wm : ion OEE EE ‘0 o*o 0 0 ao ao oo, 


Slo ZHIO Soo 


oe wo cor? 
30N).. IdM 


oo" 


oo'2t 


Y3A07 AGUINIWIC3S SNIddiO °930/01+ Ol Loy Y34SNGYL 


as 


(2) 


(pd) 


OHOh 


INI dd IO 


“930 


Ol- 


OI Ley 


Y34SNUYL 


OI LbYy 


Y34ISNBYL 


oo"mt 


f ; 
7 

‘| 

y 

eh 


|] 


ns oo ae viet “S's 


 t av t we; yey engl 


. 
t 


: ry 
7 mi v 
ay 


i 


\ =e Gtekiae SNe 


o 


— 


oe 


\ 
— 
= 


x 


70 


because it is also a function of terms used in ray 
expansion and this depends on the region in which 
rays are to be evaluated. Long period spectrum is 
very sensitive to the terms used in ray expansion 

and it is clearly also dependent on angle of 
incidence of the wave. For severely dipping inter- 
faces, short time interval or high frequencies should 
be made of use. 

Dip determination using spectral ratio method 
requires important assumptions. For even a simple 
crustal structure, the experimental transfer ratios 
for various azimuths have to be studied to determine 
the amplitude changes indicative of a possible dipping 
structure. If short period data is used to study 
neha interfaces, assumption of non-attenuating and 
isotropy of the media is necessary and velocity or the 
media and angle of incidence of the wave have to be 
known with sufficient accunacy.. “for a2 mulerulppiIng 
layered media the spectral ratio method can then be 
used to obtain a better fit between theory and experi- 


ment, if additional evidence exists for dipping layers. 
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‘Ley Representation of attenuation in the solid Earth 


The study of seismic wave attenuation in the 
solid Earth has developed recently into one of the more 
interesting research areas. In many seismic studies 
attenuation has been neglected for the sake of simplicity. 
However, a number of researches have shown that a study 
of the loss mechanism brings additional information 
about the physical state of the interior of the Earth. 
The nature of the loss mechanism is very important when 
inversion of observed data is attempted because a 
seismologist would like to interpret his findings in 
terms of a real Earth. 

It has been customary to represent a lossy 
media in terms of a delayed strain due to an applied 
stress. This can be written in the time domain as a 
convolution of strain with a memory function or in the 


frequency domain it will be 
P(w) = e(w) -m(w) ; w = 2nf (1) 


where P(w), e(w) and m(w) are the Fourier transform 

of stress, strain and a memory function respectively. 
Here the memory function, m(w), describes the nature 

of the loss mechanism such as viscosity, solid friction, 


etc. The memory must be an even function of frequency 
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in order to have energy dissipated equally over positive 
and negative frequencies. There are a large number of 
models suggested to describe the loss mechanism since 
Maxwell (1866). They are usually expressed in the time 
domain by a suitable approximation of laboratory measure- 
ments. By calculating the Fourier components of the 

time model, attenuation is also aeestield in the frequency 
domain. A large number of laboratory measurements on 
solid materials (Knopoff, 1964) showed that, for strains 


ceeroretharl won. 


, attenuation is independent of amplitude 
of strain and linear for sufficiently large frequencies 
but non-linearity occurs around zero frequency. 
According to Hooke's law, the ratio of stress to 
strain is a parameter called the modulus of elasticity. 
Then in equation (1), the memory function is nothing 
else than a complex modulus. One way of treating the 
problem is to express an elastic modulus as a complex 


function in the frequency domain. For example, in 


isotropic lossy media the two elastic parameters will 


be 
= Ap + 2 Sgn. Ww Ay 
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where A and wu are Lame constants and the sign w function 
is required so that there are real space-time loss 
mechanisms and stress components. This approach can 
be thought of as a modification of Hooke's law. As 
will be shown below, introducing complex elastic modulii 
is equivalent to having a spatially damped motion in 
frequency. In equation (2), the elastic parameters have 
been assumed to be independent of frequency. This is 
equivalent to having frequency independent velocities. 
However, Futterman (1962) calculated the Fourier trans- 
form pairs for a model in which body waves are dispersive. 
For the frequency independent case, the wave 

velocity in a lossy media can be expressed as (Bullen, 
1953) 

3 Mp +14 S2GN. M, 


VO ata a a ae : CS) 
p 


where op is the density of the media and M is a complex 
elastic modulus. As an example for P (compressional) 


waves 


M, = AR ats 2Up ’ M, = Ay ct 2uy ° (4) 


From the definition of a phase angle 
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a M for My <i Mp 4 


The inverse of the phase angle, 6, is defined as 
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"quality factor, Q". This dimensionless quantity has 
been used as a measure of attenuation in electrical 
circuit theory and many branches of physics. There 
are a number of different definitions of Q in the 
literature, however, by making use of various physical 
quantities, it can be shown that they are equivalent. 


For example, it can be defined as 
Zi 
) += AE (5) 


where AE is the energy dissipated and E is the maximum 
energy stored in the system per stress cycle. Another 
definition is that in a spatially damped motion where 
the restoring force is proportional to the amplitude 

of vibration, the natural logarithm of the amplitude 
ratio of successive oscillations (logarithmic decrement) 


is proportional to the inverse of Q 


(6) 


in ail rises definitions, an exponential type attenuation 
results. If we consider a plane wave propagating in 

the z direction in isotropic lossy media with a certain 
velocity, then at a time t = 0 it can be written as 


-iw (p/M, +i sgnwM,))*z -|w|z/2CQ =402/C 
8 =e 2 . (7) 


This represents a plane wave propagating with a velocity, 
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C, whose amplitude is damped with distance z in an 
amount determined by the frequency and the factor, 


Q. The phase velocity is given by 


2 
2 


} (8) 


1 p 
¢ M, pice ve sgnw M, 


and the spatial attenuation constant will be 


Be.) ae e % 
| aCe Nee ‘Wl Fa sgno M,) ¥ (9) 


Equation (7) is a particular solution of the differen- 
tial equation of motion of a medium possessing solid 


£rict#on” (Knopore, 1956) 


ik 3 = 1 3 — 
¢) ey ae (A+ 21) £2 Tuy eee oa We (meee, 


(10) 


where oe and Qe are quality factors for P and S waves 
respectively. 
Solid friction is one of the physical models 
for attenuation where the phase velocity and the 
quality factor Q are not a function of frequency. 
This model is strongly suggested for solid materials 
due to laboratory measurements but there is no theory 
to explain every aspect of solid friction. 
We have seen that attenuation can be sufficiently 
described by having complex velocities. Using equations 


(S) ana a ao, M, can be eliminated: 
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acre 42 i (11) 
PaO rr. | 


If Q >> 1, then approximation will be 


coi ae Giz) 


On the other hand, using same equations (8) and (9), it 
can be written 
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Then equations (12) and (14) determine the real and 


imaginary part of velocity. 


2. Attenuation properties OL Chico 


In many research works,attenuation properties 


of the crust has been treated together with attenuation 
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in the mantle or attenuation characteristics have 

been used by comparison of laboratory measurements. 
However, it is difficult to imitate physical states of 
the interior of the solid Earth in’ a laboratory. Born 
(1941) discussed the nature of solid friction and 

viscous friction and attempted to verify his discus- 
Sions by experiments using bar type sample rocks. 
Oscillatory motion has been set up at one end of the 

bar and motion is detected at the other end by a coil. 
From the damped resonant motion of the bar’, he calculated 
the logarithmic decrement using the width of the resonance 


curve. His results can be summarized as in the follow- 


ing table: 

Resonant freq. (K.Hz) noe 

Sylvan Shale D2 OU. 70 

| 6.58 5) 

Behe 73 

Hunton Limestone 10s 5-31 60 

| 2. ie ne 

Amherst Sandstone eo Ss) 

0.93 roe 


Born used the above seven dry samples and concluded that 
losses are due to solid friction but unconsolidated 


sediments near the surface may exhibit plastic behavior. 
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Therefore, there could be viscous’ losses. Similarly, 
under laboratory conditionsBirch and Bancroft (1938) 
studied forced resonant vibration of a long column 

of Quincy granite. They used longitudinal, flexural 

and torsional waves to measure dependency of velocity 

and attenuation on frequency. Their result showed a 

Q value between 100-200 over a frequency range 140-1600 
Hz. However, the same authors (1938b) measured a Q 

value of 1200 for torsional vibrations of a small granite 
Foa™at avrrequency about FLO" KiHz“under*a pressure of 
4000 atm and a Q value of 180 for a pressure of 200 atm 
at the same frequency. Knopoff and Porter (1963) deter- 
mined Q = 80 for Westerly granite from the Rayleigh wave 
analysis in the frequency range of 100-400 K.Hz. There 
are a number of research work done in measurement of 
attenuation in situ. Among these, McDonald et a1.(1953) 
have done a detailed study of attenuation in the Pierre 
Shale near Limon, Colorado. They measured attenuation 

of vertically travelling compressional waves using 
signals recorded by five detectors positioned in the 

bore holes at various depths in the formation, which was 
about 4000 ft. thick. An empirical frequency decay law is 
stated using the frequency spectra of recorded signals 

in the. frequency range 50-500 Hz They concluded 

that attenuation is a function of the first power of fre- 


quency contradicting the work of Ricker (1953) at the same 
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site, who found an attenuation relation which is a 
function of the square of frequency. The McDonald et al. 
result showed a Q value of about 10-23 for Pierre 

shale and no dispersion observed. However, they have 
noted that dispersion may be important in the weathered 
zone. Collins and Lee (1956) used a limited number of 
strain recording taken over a small distance to deter- 
Mine attenuation parameters in sandstone. Their 
results gave a Q value of 7 over the frequency range 100- 
900 Hz for compressional waves. However, they could 
not make any conclusionsdue to insufficient data. 

Many laboratory measurements for attenuation 
have been performed at large frequencies. One obtains 
a high resolution in the data and so determines the 
dispersion relation with frequency. It is well 
agreed among a number of authors that attenuation 
is linear for these high frequencies. However, 
the seismological spectrum is centered at about 1 Hz 
which is well below the range of laboratory measure- 
ments. On the other hand, it is well known from expe- 
rimental seismology that waves with a frequency near 
50 Hz create great problems in interpretation since 
at these high frequencies elastic waves are very 
sensitive to inhomogeneity and anisotropy of the crust 
and scattering becomes important. The first low fre- 


guency attenuation experiment has been performed by 


sLxerg tot Ecos 
eved yes ,yevewol a bee Te 
bersdtsaw ods gut Rae em vei notixegeth tnt 


gis 
to secu Sothmtl .s bean (age), ed bas. ami L09,, +6 
riedeb of soastetS~ i sent & Teve netss | _ ReiLbtone, a. 


afeit .enotebuse at a 

~O01 epacx, Youenperl odd aeve /¥ ae autor 0 * evap ¢ 

hives yet \.ovewol .pevew tendienongmne atene 

-etebh tapioitivent ad UE BA ERLE FE 

noikeude3ts 101 ei nepisivesse: protiuedsl seni 7 

enisddo ef60 .eeionevpest eptset 26 bemxotzeg 2 
oft eoninreteh og bap stab eft nk aided 

‘iow ei #1 .youeupett dtiw in: 

noiseunesis ted exodiup, ‘ke sedmuie, a ‘proms. b | 

‘tovewoH .eatoaoupest pid seatt xo? incisal 2 ' 

sill sods 46 botedads et jurxsveqa\isntpoloma be | 


<0 


B af i 
- " : 


1 


Kp 


~sisesam yrotsiodel 20, epnsx odd esta mre at asia | 
-agne Moz? awonl flow ab dk «bes asdto eds, 10 sateen 
2a yonpupaxt & Hatkw Bovey steris ehosommtos tegnamis |.) 
ecaiche rolissoxgyodai at emmlderg seorp saeers AME | 
YXOV gis eovew oitashe #8to0 
tewr edt lo yqoxtesian. - 3: 
“932 wol jexli sdT tkoge 


80 


Peselnick et al. (1961). They have attempted to measure 
internal friction in shear modulus of Solenhofen lime- 

stone over a frequency range of 4-10' Hz. Using torsional 
waves for a MEE Cone sample whose density is 2.67 gr/cc 


and shear velocity is 3.1 km/sec, their results showed 


higher Q values: 


Solenhofen Seog Hz 850 


Later, Donato et al. (1962) performed attenuation 
measurements in the laboratory using a wide range of fre- 
quencies including very large periods. For low frequencies, 
the specimen was made to vibrate as a beam clamped at both 
ends and then compared with a standard steel beam. 

Young's modulus and attenuation have been measured. 

They found that while Young's modulus remains constant 
within 5% for all frequencies, the Q varies at low 
frequencies, but it is constant for high frequencies. 
Variation of Q at low frequencies is noted to increase 

by a factor of two as compared to the value at high 


frequencies. Their results are as follows: 
ni Q 
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£ Q 


Dolerite Sie 4 HZ 360 
2-40 Hz 180+18 

Hard Chalk <2. Hz 2:20 
11-66 K.Hz 70 + 8 


Tullos (1969) has used transfer function analysis to 
measure attenuation in Gulf Coast sediments. Analysis 
has been carried out using vertically travelling waves 
in a sedimentary layer of 1000 ft. thick, recorded in 
a frequency range of 50-400 Hz by array of seismometers. 
Q values have been found to be high for these sediments 


for which the velgcity structure is well known: 


d (£t) Qp 
Loam-sand-clay 20 5 
Clay-sand 75-100 420 
Sandy clay 100-500 165 
Clay sand 500=F000 yh, 


There is a large amount of research on the 
attenuation of the Earth's interior. In those works 
where researchers treated the crust as a part of the 
mantle crustal attenuation is found to be severe. 

There iS a small number of published works dealing 
with attenuation only in the crust. Press (1964) deter- 


mined Q values for the crust from nuclear explosions 
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in Nevada. He made use of channel waves and assumed 
that seismic wave energy decreases with the first power 
of distance. His average for Q from Ly waves is Q = 

450 and Q = 260 from Py phases. He interpreted low 

Q values for at waves not as a property of crustal 
structure but as an energy loss due to mode conversions. 
Fedetov (1963) studied attenuation of shear waves 

whose frequency is 0.2-0.5 Hz propagating in the crust 
near southern Kurile island. He found a Q value with 

a large uncertainty: Q, = 400+200. Anderson et al. (1965) 
estimated a Q value of 1000 for the crust using surface 
waves with periods 50-300 seconds. O'Brien (1968) 
estimated a Q value of 1000 for compressional waves 
using Early Rise refraction data near Superior Lake. 
Sumner (1967) made use of earthquake shocks located 

at a depth of 20-180 km, to study attenuation of P waves 
in southern Peru. He calculated a Q value of more than 
1000 for the crust and upper mantle. Clowes et al. (1970) 
determined crustal attenuation from deep reflections. 

His results showed a Q value of 200 for sedimentary 
layers and 1500 for the lower crust. Hill (1971) used 
published head wave amplitude data to determine Q values 
in central north America. His results from amplitude 
distance curves of Pg and P* arrivals, gives a Q value grea- 
ter than 1000 for the crystalline crust and intermediate 


layer in the west coast normal heat flow provinces and 
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Q values are smaller for high heat flow provinces. 
He attributed low Q values to the scattering of head 
waves, temperature dependent Q values and negative 
velocity gradients in the upper part of crystalline 
crust associated with high geothermal gradient. In 
his work, he used a crustal model with 2 km low 
velocity sedimentary layer (Vp = 3.0 km/sec), overlying 
crystalline basement (Vp = 6.0 km/sec). 
| Brian (1973) utilized 4-40 sec period surface 
waves to determine attenuation in central north America. 
He estimated an average Q value for shear waves of 
about (75-300) for upper 15-20 km portion of the crust, 
Q = 2000 for deeper part of the crust. He interpreted 
this rapid jump in Q values as a possible temperature 
increase and change in chemical constitution of the 
crust at intermediate depth. However, he was not able 
to. predict a,Q value for P waves. 

In a tectonically complex area, O'Neill and 
Healy (1973) reported very low Qp values along the San 
Andreas are region. For the uppermost part of the 
crust they have estimated Q values of 35-100. In 
some region Kurita (1973) has elucidated Qq, values 
about 90-160 for the upper 6 km of the crust and for 


sedimentary rocks it is as low as 35-40. 
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te epectrel ratios for dissipative crustal models and 
application to earthquake data 


In one of the attenuation models, the memory 
function in equation (1) is expressed as a logarithmic 
function in time which fits a linear behavior of atten- 
uation over a large range of frequencies. In this model 
phase velocity and quality factors are functions of 
frequency (Futterman, 1962). However, a cutoff fre- 
quency is introduced to approximate laboratory 
measurement. Above the cutoff the Q is linear with 
frequency, while below the cutoff frequency, Q is assumed 
to be infinite. Wuenschel (1965) and Jensen (1970) 
applied this model to real materials and spectral ratios. 
In this thesis, frequency independent Q value and phase 
velocities are used. Attenuation is introduced in the 
Haskell-Thomson formulation by making use of complex 
velocities. Attenuation is allowed to be independent 
of the sign of frequency and distance travelled by the 
wave in the z direction. The proper sign is chosen for 
the exponent to account for the correct direction of 
transmission for a wave incident at the bottom of 
layered media. In Figure 1, the transfer ratios for 
one of the isotropic Alberta models are shown for 
different Q values. In this figure, the solid line 
represents a non-attenuating model. The attenuating 


models are only shown by their envelopes on the graphs. 
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Figure IV.1 


Comparison of transfer ratios for attenuating 


models with the non-attenuating one. 
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Very small values of Q cause the spectral level to 
increase with increasing frequency. For example, for 

Q = 50 significant differences from the non-attenuating 
model “occur at about 2 Hz, 2.60 /Hz and 3:40 Hz. 

Among these, the one at 2.60 Hz is the most indicative 
for attenuation. This peak has an 800% larger amplitude 
than the second larger peak at a frequency of 3.40 Hz. 
Models for Q = 500 and Q = 2000 have very little dif- 
ferences in amplitude from a non-attenuating model at 
all frequencies. Letting Q increase with depth in 

the crustal section (Figure 2a) creates no distinct 
differences from a non-attenuating model. All models 
have negligible differences in amplitude for frequencies 
tess than “0.90 Hz for 35 degrees of angle of ancidence. 
For angles of incidence of 5°, all models are almost the 
same at frequencies less than 0.3 Hz. Among these 
calculated models, dissipative ones have shown no 
frequency shift for a sampling frequency of 0.0122 Hz 
with respect to non-dissipative models for long periods. 
But major peaks which are indicative of attenuation at 
higher frequencies are displaced as much as two samples. 
Here, it is assumed P and S velocities have the same 
quality factors. However, if we were to use a different 
quality factor for S waves following Knopoff, 1971, it 


would be 
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6 a T 2c (14) 
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Figure IV.2 


A transfer ratio for attenuating model 


and its time synthesis. 


87a 


TRANSFER RATIO 


Pp 
200, 2000, 2000, ©, . ee daveb 


6.40 


VERTICAL i=35° Q=200, 2000, 2000, %, 


AMPLITUDE 
3.20 


RADIAL 


13.00 


12.00 
SECONOS 
RADIAL : 


oo 


GV G ro) 


iy raat ae) 


Ba ee 


HLT TAME 1 3 
Pal te & (at s. ee consioobe ey ve 


or aa “ea ae 


en a Ee eee Teed ee ee. as ~~ — oes 
a i TPRART AE Pe 
ie; 


: F * 
a) ia ees 
‘ 


: if 
,. 


dodnae i tis ioe 


a + stipe ii 


az \ 
~ 
' 7 
y ! 
¥ 1 Diy 
. | 
if 
, ry 
i . } 
i MM 5) 1 
ery t 
é ® t 
ae ob a 
ra p - ; Tie 
2 ) ; 
y ies y i eo . 1 pal n 
i ‘A _— 4 iS 


88 


then more distinct amplitude changes will occur in 
the transfer ratios. This is due to exponential 
Eactor: 

neers a 
where Vo< Vp and Qn? Qa. Then the exponent will be 
larger than the similar exponent for P waves for the same 
frequency and distance z. Therefore, the horizontal 
component will be damped faster than the vertical com- 
ponent for all frequencies. 

A synthetic seismogram Late studied for an 
attenuating model by taking a Fourier inversion of 
the transfer function, provided that the true phase 
is known at each frequency. In Figure 2b, a synthetic 
seismogram is determined for an Alberta model including 
attenuation. There are significant amplitude differences 
from the non-attenuating model. The time synthesis was 
obtained by a convolution of a pulse with the crustal 
transfer function calculated using Haskell-Thomson 
matrix formulation. The source pulse is a sine func- 
tion which has a decaying envelope in time and the 
frequency of the pulse is centered at 10 Hz. More 
detailed calculations of theoretical seismograms will 
be discussed in another chapter. What is clear from 
Figure 2b is that the amplitudes of different phases 


are decreased with respect to the non-attenuating model. 
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However, the relative amplitude decrease among phases 
is nearly constant. This makes determination of a Q 
factor, by making use of the amplitude ratio of two 
or more phases in time domain using experimental 
seismograms, very difficult. When very low values 

of Q are used, the shift of many phases is about 0.3 
sec in the time domain. This is because large phase 
shifts are introduced in the frequency domain when Q 
values are small. Large phase variation due to severe 
attenuation might make Fourier inversion unsuccessful 
when time synthesis of a model is required. 

When determination of attenuation is desired 
for a crustal structure using transfer ratios, exact 
amplitude correlation between observed and theoretical 
responses is necessary. This requires a well calibrated 
recording system. On the other hand, for a given 
crustal structure, if there is a sedimentary low velo- 
city layer, velocity and thickness for the layer must 
be known accurately and velocity contrast between layers 
has to be reasonably known in order to interpret large 
amplitude variation in the transfer ratios. Never- 
theless, modelling of large amplitude variationsas in 
Figure 1 for the case of Q = 50, by using velocity 


contrast or various Poisson's ratio is difficult since 
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these quantities generally have an upper and lower 
limit for a@irealistic:-Barth's crust. On the’other 
hand, due to finite experimental data, a calculated 
spectra using real data is convolved by a kind of window. 
Generally, this window is not a simple one in order to 
increase reliability of the calculated spectrum. For a 
better comparison of theoretical and experimental spec- 
trum amplitudes, similar windows must be applied to both 
spectra since the window is a significant factor for am- 
plitudes as well. A spectral window less than 50 seconds, 
eliminates completely large spikes as in the case of 
Q = 50 in Figure 1 and makes the transfer ratios indis- 
tinguishable for various Q values. However, for spectral 
windows of more than 50 sec length, spectral amplitudes 
are independent of window length and only in these cases 
more convincing comparisons of amplitudes between theor- 
etical and experimental spectral ratios are possible. 
Another significant factor affecting the ampli- 
tudes of experimental spectral ratios is additional 
phases such as near source or core reflections. These 
additional phases create multiple pulse problem ina 
combined layered response. Even when these additional 
phases arrive with same velocity,spectral ratio amplitudes 
are found to increase significantly. If deep earthquakes 
are used for the study of attenuation, then additional 


phases can be neglected in the analysis since these 


‘alli bab teqai ne é 
tetgo ot? 2 aves : vk 
bets tooten & . s4a6 Lahets : rhe ‘ >  basd 
Wohitiw Io bart &)\ yd: poviovact ‘ak "ha da eakaw “exd8 Hoe > 
0? 1860 ai one efgate’ 6s Soa di wobatw aide aye wer ; 
5s 10% jmbttoeda besaluatne edt Xo lqottbdalion | a ee 
-pede Latusmixegxe bie ten tteroods to ‘foabisqmon a 


‘itod of Beitlaaa sd seum awobaiw sekbatte ‘Bebutita 


Jed es 
~ms i692 Yoshet snesitingid a 2 wolbatie ent ‘esnke 
,ebnover OF nett acol wobdiw ferivega A ‘few as 
et At 

to sano ofa ni 2h asdiqa spraet eisdelgnion ‘eats 

” a 
; ~ Au 
-2ibat eok3 étensxd odd Seinwst bas i abet t ah 4 he 
lextoeq@a tot \sevewoll , setiisv 0 ari and 103 otdateko A 
: Net <a 

esbutiignuse L[sxtosqa ,dgpaet sen 02 ued 9 Sxom hud wobs 
e08s> seodt ni-ylac bas depael wobrtw 30 


-xeadis tsewied. sebytifgqms Lo anoadtedeea: R 
Cm sin 2ottsx Lextooga istnembxegne Sl 

~Ligqne, siz pabtoetie 4ogosn? sacsthiaghe sedzona 

isnoithibts ei sokatnt Iszassge Istnembteqxe — . " 

seol?T .enoltoelie:s S105 HO sbawea ‘tent 26 foue eh 
5 a Metairie, selog aligétton baasindetd ney, kasi 

* Ianoisibba auieats Role moved ‘ee yet £ Dents es 
sobtiigns obtes” Low sroage .yaiteatiay os — whos ae A A 
sanliateniitl ideal +x nchsiibothamge 8 ios ‘0 {os | + bayot exe ee ey “ 
fsaoitibbs aed mod teurtedts > hase Pree sot beau ox8 | Bey 


)) Pabdd ovate etayteas Seth alata od abs esesdg Me Ok 
: a | ee Z oc oO “, 3 *) Ae RAL, 


(et i ' " 7 : - Aes 
ee : Vs is wear 


_ 


Pi ah 
a f iy 
Nive care . 


ei ne 
Le fe 


ir. be Layee. Cb ar vg vrs ms "a . anc) 
ON) Ge Me hime TR 


A Ae a ee 


» tp Cae 


Figure IV.3 


An experimental transfer rdtios for a large 


window length. 
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phases arrive out of the window length for earthquakes 
which occurred at far distances and large depth. 

Effect of additional phases and earthquake magnitude 

on the amplitude of spectral ratios will be discussed 
extensively in another section of this thesis. During 
the spectral analysis of 1970 earthquake data recorded 
by the VASA stations in central Alberta, an attempt 

was made to determine attenuation. Several earthquakes 
had no additional phases incident within 50-70 seconds 
of time. One of the results recorded at Pine Lake is 
shown in Figure 3. This seismogram has no predicted 
additional phases up to 102 seconds when the pP 
reflection arrives. The core Perieccicns are inci- 
dent at the same time as the main P arrival but this 
has no effect on spectral ratios (this will be shown in 
another section). Spectra of vertical and radial motion have 
been calculated by inversion of the smoothed auto-correlation 
function. Smoothing has been achieved by a 70 second 
Parzen window. The earthquake was located at 77.6° 
distance away and the predicted angle of incidence at 
the bottoem-of-cerust-1s about’ 25-30°, corresponding; to 
a phase velocity about 19.25 to 16.16 km/sec. This 
spectral ratio was comparable to the theoretical one 


for a model shown in Figure 2a. Distinct features of 


the model for different attenuations are almost kept 
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unchanged. If we study the experimental ratios shown 

in Figure 3, there are no features resembling the model 

with Q=50 but agree fairly well with the one with Q>500. 

For this interpretation,assumption of constant phase velocity 
and the well known velocity variation with depth as 

shown in the model is necessary. 

Spectral ratio models do not discriminate 
against various attenuating models for a source 
distance less than 16°. However, aie ferences between 
non-attenuating and severely attenuating models are 
Still significant for small angles of incidence. 

From the study of peak positions in 
experimental and theoretical spectral ratios, it 
is clear that there are significant differences. 

Perhaps it would be possible to say more about atten- 
uation when a better fitting model is found at higher 
frequencies and if gains of recording instrument and 
other information are available, absolute match of 


amplitudes may be attempted. 


4., Discussions and conclusions 


From the review of published research, a number 
of conclusions can be reached. Laboratory measurement 


of attenuation cannot be accepted as a representative 


value for the real Earth. The consensus of opinionis that 
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for high frequencies ,attenuation factor is linear. However, 
nonlinear attenuation exists around zero frequency and 
the exact cutoff frequency is not certain. Unfortunately 
the range of our seismic spectrum lies in this nonlinear 
region and there is a larger scatter of data in this band 
of the spectrum. Mechanism of loss in solid materials 
and tne Earth is still an unresolved problem. Today 
one of the two principal theories proposed to explain 
anelasticity is solid friction, where Q and phase 
velocities are not a function of frequency. (Knopoff, 
1956), the other is causal dispersive attenuation 
(Futterman, 1962). There are a large number of experi- 
mental data supporting non-dispersive attenuation in 
solids. However, Wuenschel (1965) showed that very 
small dispersion can significantly affect attenuation 
characteristics. Both theories have a model for the 
nonlinear attenuation region, but no firm experimental 
evidence exists to support any of these. 

The nature of the loss mechanism in the real 
Earth can be divided in two parts: 

a) Losses due to rock formations in natural 

states. 

b) Losses due. to other causes, 

In this classification, part b)covers a wide range of 


causes such as reflection, transmission, scattering, 
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etc. In general, this is a representation of the 

real Earth in terms of wave propagation. Part a) 
represents attenuation, for which measurements 

have been attempted, that is stress. relaxation 

due to the atomic or macroscopic structure of 

rocks. In an experimental attempt on the real Earth, 
these two causes generally are inseparable. For :a 
correct measurement of the Q factor, the Earth has to be 
modelled to a sufficient degree of approximation. Channel 
waves and head waves, due to their propagation patterns, 
are sensitive to the lateral inhomogeneities in the 
Earth. Attenuation measurements using these waves, 
generaliy, show large variations even for the mantle 
region of the Earth. 

Q values increase with increasing pressure and 
temperature but after a certain temperature, thermo- 
friction becomes important (Knopoff, 1971). It is most 
likely that the attenuation mechanism becomes more com- 
plex in the deeper paLt (OF (they cyusc. 

A solid friction loss mechanism predicts a lower 
Q values for S waves. This is due to the argument that 
permanent creep is almost never observed in a deforma- 
tion experiment in which the sample is compressed 
hydrostatically but not fractured (Knopoff, 1971). 


This argument is against the other theories of loss 


mechanism (Walsh, 1969). There are a large number 
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of measurements which showed iargee Q than Qp (White, 
1965). On the other hand, from surface wave study, | 
Anderson et al. (1965) observed that Q values for P 
waves are 2.4 and 2.6 times larger than Q values for 
S waves. 

In view of the experimental data mentioned 
in section 2, it will be fair to assume Qp = 200-250 
for a 2 km thick sedimentary layer composed of a series 
of shale, sandstone, dolomite, limestone, etc. and 
Qp = 1000-1500 for a granite basement and lower crust 
in Alberta. During the analysis of 1970 earthquake 
data recorded by the VASA stations, experimental spec- 
tral ratios have been compared with various dissipative 
theoretical transfer ratios. This comparison suggests 
an average value of Q, 2 500. However, the results are 
not conclusive due to an insufficient amount of reliable 
data. More reliable Q values as a function of depth 
cannot be determined at this time, using the spectral 


ratio method. 
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V. Crustal Structure in Central Alberta 


from the Spectral Ratio Method 


dre Experimental data 


In this thesis, an attempt is made to resolve 
crustal structure in Alberta using about 50 earthquakes 
recorded by various seismic stations (Figure 1). The 
data were recorded on a magnetic tape digitally with a 
sampling rate of 12.5 samples per second. Spectra 
were calculated using the method mentioned in section I. 
Experimental spectra have been compared with the two 
dimensional theoretical ones. Azimuthal deviations have 
been measured at every station to determine the effect 
of wave scattering in a three dimensional case. Station 
residuals have been used in the interpretation of spec- 


tral. ratios. 


2. EDM station 


Event 73 recorded in 1970 has no additional 
phases predicted within the window length. This event 
has a prominant peak at 0.1 Hz (Figure 2). Other 
events have a peak at about the same frequency but 
small differences in peak position exist. Events 71, 
92, 93 showed anomalous shifts in peak position at 


lower frequencies indicating a thicker crust. However, 
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Figure V.1 


Location of array of ‘Stations. 
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Figure V.2 


Long period spectral ratios for EDM station. 
Frequency of expected holes in the spectral 


ratios has been indicated by the arrows. 
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all these events suffer distortion due to additional 
phases such as pP and PcP. These additional phases 
produce holes in the spectra, which distort the ratios, 
but the large observed frequency shift of 10 mHz 
cannot be attributable only to the additional phases. 
Short periods of spectra were found to be more 
sensitive to the noise. By noise we mean signal gen- 
erated noise such as anomalous conversions in crustal 
layers or surface wave type noise. A large peak at 
about 0.5 Hz was very well duplicated on all events 
(Figure 3). A second peak at 1.60 Hz was also well 
established but it has a lower amplitude. A third 
peak at about 2.5 Hz was not seen from all the events. 
However, according to event 73, it must certainly exist. 
Averages of all spectra were used to display in 
Figure 3 the most probable peaks and their positions. 
In Alberta a large number of crustal investi- 
gations have been carried out in the past. However, 
most of these studies have been done in the southern part 
of the province. In earlier spectral ratio studies 
(Alpaslan, 1968; Sprenke, 1972), a southern Alberta 
crustal model (ALTACRT-2) has been used as a starting 
model. This model did not fit the observational 
spectral ratios as shown in section III. A reflection 
study (Ganley and Cumming, 1974) near EDM has brought 


some clarification to the crustal structure in this 
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Figure V.3 


Short period combined spectral ratios for EDM station. 
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region. Although this work is not complete, clear 
reflections have been identified at 3.5, 7.5 and 
10.7 seconds. According to this work a model 
has been constructed for this area (Table ITI) 
which has a total thickness of 36 km. A low velocity 
layer of 2.3 km has also been included. In Figure 4, 
possible models have been displayed. A low velocity 
layer did not have a pronounced effect in this band 
of spectrum. Removal of the sedimentary section does 
not influence peak position but significantly affects 
the shape of the spectrum. By removing a 4 km thick 
intermediate layer (model 3), the peak position is 
displaced by about 4 mHz towards the lower frequencies. 
This is an expected phenomenon because the high velocity 
layer is displaced by a lower velocity layer which 
decreases the travel time. In this frequency band, 
layers with velocities of 6.4 and 6.5 km/sec are not 
distinguishable. In Figure 2, the experimental spectral 
ratios were compared with model 2. It is clear that 
the peak of event 73 has been slightly displaced to 
the lower frequencies. If we keep the crustal thickness 
constant at 36 km the experimental spectrum favors a 
simple layered crustal model. 

Then, is there really an intermediate layer with 
a high velocity? Same model can be produced if we have 


thinner crust with a low velocity or thicker crust with 
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Figure V.4 


Long period spectral ratios from models given in 
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Table II. 
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a high velocity.as in Figure 7. Residuals at this 
station show an early arrival at all quadrants. 
Therefore, a thicker crust with a high velocity zone 
is possible but the amount of the residuals does not 
suggest any large changes in thickness-velocity rela- 
tions. 

The sedimentary layers of the model control 
the shape of the spectrum at higher frequencies. 
The deeper layers of the model have an influence on 
the amplitude of the spectral ratios. A two layer 
sedimentary structure has been found to be satisfactory 
to produce the basic pattern of interference which is 
observed (Figure 5). In this figure, event 71 has a 
peak of very large amplitude at about 2.5 Hz. This 
is due to the presence of additional phases. Dipping 
models produce no significant shift in peak position 
with respect to the flat model but only the amplitude 
of the peaks changes. A very thin low velocity layer 
does not have significant effect on the spectral ratio 
at this frequency. Therefore, this band of the spec- 
trum is not sufficient to resolve a thin LVL as suggested 
by Somerville and Ellis (1972). This conclusion perhaps 
is not Surprising that 'ayvthin ‘layer iike ‘this LVL ‘can 
be resolved for wavelengths less than layer thickness, 
that is a frequency band larger than 2.5 Hz should be 


used. Small peaks in Figure 5 do not compare with the 
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model as well as large peaks; a better sedimentary 
model can be found when larger amount of reliable 
spectra are available. In many ways amplitudes of 
the peaks can yield more information, but for this 
kind of comparison all information regarding instru- 
mentation must be known. Inexact gains of experimental 
spectrum component have made amplitude comparison 
doubtful. Record character for this station was 
poor (Figure 6) in the sense that additional phases 
change the appearance of record as well as the frequency 
content of it. ArrivalS with a larger amplitude than 
main P arrival is a general character of the records. 
As can be seen, large transverse motion for 
all events is a common factor. This suggests that a 
two dimensional model is too simple for this area. 

If the crust were a homogeneous media transverse 
Motion would be zero. However, inhomogeneities will 
produce transversely polarized shear (SH) motion which 
should have a minimum in same direction. Starting 
from this assumption horizontal energy ratios were 


calculated (Kurita, 1972): 
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Figure V.6 


A few examples of records for EDM station. 
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Figure V.7 


A graph of azimuthal deviation versus azimuth 
and station residuals. For this diagram‘’all 


available events were used. For residuals J-B 


tables were used. 
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where ER and E, are the radial and transverse amplitude 
spectra respectively. Using the known station-source 
azimuth as a starting angle, data were rotated to 

obtain the minimum on value which can be considered 

as a true wave approach. The sum in equation (1) was 
performed from 0 Hz to 4 Hz to obtain the minimum 

energy for the transverse component. The azimuthal 
deviation was measured for all available earthquake 

data. As seen in Figure 7, a large amount of variation 
exists. However, one should remember that this method 
includes possible misplacement of seismometers and misloca- 


tion of the epicentre, Therefore, the pattern of azimu- 


thal deviation is more important than the magnitude. 


3. “=PIN statwvon 


For this station, data free from predicted 
additional phases were found only for the SE sector. 
Events 48 and 68 yield a peak at about the same 
frequencies but 10 mHz difference exists in peak 
position. Event 73 has very low amplitude and did 
not produce a peak as other events. These differences 
cannot be explained using the theory discussed in this 
thesis. However, event 48 has a maximum amplitude 
and when combined spectra are joined (Figure 8), SE 
group produced a peak which has the same location and 


shape as in EDM station. Model 2 of EDM produces a 
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Figure V.8 


Combined long period spectral ratios for PIN 
station. Event 83 has no resolution at low 


frequencies due to lower magnitude, therefore 


not included in the averages. 
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very good fit to the observed combined spectra (Figure 
8). Events belonging to the NW sector produced very poor 
results, the shape of spectra being distorted and shifted 
to lower frequencies. All events for the NW sector 
contain additional phrases (Figure 8) and no dependable 
event has been found. 

The short period band of the spectrum (Figure 9) 
has been found to be unstable. Peaks in the spectrum 
have shown variations from spectra to spectra. Peak 
at about 0.8 Hz has fairly well been produced by 
different events and tne double peak at about 1.5 Hz is 
also a common feature of the spectrum. Comparison with 
EDM short period data shows that at this station the 
sedimentary layer is much thicker and this layer is much 
more complex due to a factor affecting the long period 
of the spectrum. Close peak position of the experimen- 
tal spectrum suggests a thicker sedimentary layer. 
However, modelling of this spectrum using the theory 
mentioned in this thesis is aericone Plane layering 
is the most likely one of the assumptions which has 
been violated. Rapid variation in the top 10-15 km 
region of the crustal sequence is apparent in all 
azimuths. Events from the SE sectors have relatively 
early arrival times and large azimuthal deviations for 


this quadrant (Figure 10). The anomalous variation of 
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Figure V.9 


A few short period spectral ratios for PIN 


station and comparison with same models. 
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Figure V.10 


Azimuthal deviation as a function of azimuth 


and residuals for PIN station. 
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the spectrum in this quadrant should be related to the 
scattering which seems to affect the wavelengths of 


about 30-60 km. 


4. RM stations (RM1 and RM2) 


Long period spectral ratios at these stations 
have been characterized by a peak at centered around 
O.18 Hz. In Figure 12, event 73 is cléar’ from addi- 
tional phases and shows a well resolved characteristic 
peak and a small peak at about 0.25 Hz. A large number 
of events showed a peak around 0.18 Hz (Figure 12-a,b) 
but not a second peak. Therefore, we will consider this 
peak as a real one. Event 73 exhibits a very low reso- 
lution at 0.09-0.1 Hz and shows a very small peak. 
However, other events generally show a larger peak at 
these frequencies, but all events other than event 73 
have additional phases. On the other hand, event 73 
has a lower magnitude which may cause lower resolution 
at these periods. Therefore, a small peak centered 
at..about. Uvl. Hz may alsSo be accepted as a ’real one 
(Faqure. T3=—c)-. As can be seen in Figure 13-b, 
major peak positions can vary more than 20 mHz for 
different earthquakes. Events 73 and 86 produced a 
peak at similar frequencies but event 84 has the same 
peak shifted to the higher frequencies. Events 84 


and 86 have the same characteristics (Table I) 
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according to the earthquake data; additional phases 
are incident at the same time. For these earthquakes 
source mechanism has not been studied but record 
characters are similar; event 86 has 1 km deeper 
depth_of focus. However, the large standard deviation 
in depth determination does not yield a solid compari- 
son. Distorted long period spectral ratios due to one 
or more additional phases have generally lower amplitude 
and when combined spectral ratios were formed, general 
characteristics of the spectra were retained. 

The crustal structure of regions close to the 
Rocky mountains is believed to be very complex. From 
the spectral analysis it can be seen that these stations 
produced completely different spectral ratios than at 
all other stations. In an attempt to model these 
peculiar spectral ratios using available theoretical 
tools, a different starting model was used. This model 
produces only one major peak for long periods. Parameters 
of the model were acd to obtain a spectrum similar 
to the SUpesinencal one that is a small peak at 0.1 Hz 
anda larger one at 0.186°HZ. First? the’ thacknese¢s of -the 
layers were varied in turn until spectral ratios did 
not improve to give a better fit between experiment and 
theory. Next BP and S wave velocities were changed to 
improve the comparison between observed and theoretical 


Spectral ratios. This order ol variation in. parameter 
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is hot justified. +However, such an order: of variation 
in crustal parameters is necessary for the inversion 

of data using an empirical method. Model RM1 shown in 
Table III produced two peaks similar to those observed 
in the spectral ratios. Howevér, the first peak at 

0.1 Hz has a larger amplitude than the second peak. 

This model has a crustal thickness of 35 km. The simple 
layered crust as proposed for other stations does not 
produce a spectral ratio which has two peaks as in 

this experimental case. The amplitude of the second 
peak can be increased in two ways: first by lowering 
the P wave velocities only as in model RM2, which 
Gorresponds to a Poisson's ratio of 0.2 forithe whole 
Grust), #secondty, by using=a thickér sedimentary section 
(model RM4). It is well known that the sedimentary 
section-ofthis-region-is-very “thick sin fact, a 2.8 km 
thick sedimentary section increases the amplitude of 

the second peak significantly to give a better fit 
between theory and experiment. In Figuresl14 and 15, 
various models have been compared with experimental 
ratios. Events from the NW sector have larger major 
peak amplitudes than those from the SE sector. A 
dipping model like RM3 increases the amplitudes of both 
peaks. A model with extreme dipping interfaces has 


not been calculated because the spectral amplitudes 
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Figure V.1l 


Event 73 recorded by RM2 station and its long 


period spectral ratio using Hemming window. 
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Figure V.12 


Some of long period spectral ratios for RM stations 
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and their combined versions with respect to azimuth. 
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Figure V.13 


Two events from SE sector and their comparison 


with models in Table III. 
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Figure V.14 


Three long period spectral fatios and their 
comparison with RM models. Amplitudes in 


these graphs are true amplitudes. 
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Figure, V.i5 


Short period spectral ratios for RM stations and 


their classification with respect to azimuth. 
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FiguretyV.16 


Azimuthal deviations versus azimuth and station 
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residuals. 
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are contaminated by energy from additional phases. 
Amplitude differences for these two sectors can be 
explained in terms of differences in the thickness 
of the upper layers. 

The spectral ratios at higher frequencies for 
those stations show closely spaced peaks similar to 
the PIN and TRO stations indicating the presence of 
a low velocity layer or a thick sequence of layers 
in the crust. Most of the energy has been concentrated 
in a small frequency band. Peaks are broad and the 
splitting of peak occurs frequently (Figure 15). 
Modelling of these short period spectral ratios are 
difficult using a plane layered model. However, a 
sedimentary section 3-3.5 km thick is indicated by 
the model experiment. 

A number of seismic measurements has been 
carried out across the Rocky mountains, but most of 
them cover southern B.C. and Alberta. White and 
Savage (1965) estimated that there is a 31 km thick 
simple crustal section under the mountains. Gravity 
anomalies for a strip across the southern part of the 
Canadian cordillera suggest (Stacey, 1973) a thicker 
crust if this corresponds to the depth of compensation. 
As far as we know there are no deep seismic measure- 
ments around the RM stations. This region shows a 


high magnetic anomaly and a large positive Z residuals 


BS 


.enedg enol 2tbbm: mo : ai : 


Ri 1 hy 
Ai) ) Ceres 
Jt ee A 
° A 28 
; 
=r 


ad neo exez092 ows ae 


aeansok ag ‘otly fh 


coms 


toh aeinasupest xedpid ts beads: Ksadoege ont 
ot *eLimre siteeg beosea yigaols wore enolsaie « 
te soneeaxg odit ‘erriteg RBM enchiete oan bas mre 8 
exayet to eonespes ALA B 20 a etootay 
hactsettisomco med asd yersns ait to Jectt Jtegap 
sit brs beoxd ‘ot eased atm Ha 
.(@L exupst) yitnenopent earooe Asey 2 stale, i, 

sxg Boies Lautoeage boixeg trorte seeds ~26 a eet * o : 
2 ,revewoH «lebom bexeyel onska o aaken. | sivoia = < 
yd bersotbmt al Aorta wel e. E-£ noivea ‘edn — “a 


es 

nged esi ainenexuesom oimeiee 20. sxedmua A | 

to teom tud (aniadnmont ysivosl ae BBOTOR: ‘ape ioe Soa 
brie eek setgedlA baw .D, a edits seve may” 

Apids mA L6 sal attend todd beteuites (eaen) 

Yiiver. .astiedison ort <a agtiose Ledeuto 

ait to txsq atedsiog ud aeorye qiate #8 to coitbmone. 

tavtatds «2 (EVEL von de) saopgue prallibzoo as 
noLsaensgmop to ddqeb acy oF) ebnogesss0> aids 34 sem 
-Ssxvesem oimekee. geeh of ots giads wood are as. apt A 


ss 


& ewore nolper ade <200L9S2a ME. att owes sna: i Tae oe: 
| | 


afaubies: § svidisog aeuat & Bre Saas olvenpan ae a ares 
et coe ; 


V29 


on the geomagnetic components (Haines et al., 1971; 
Berry et-al., 1971). Therefore, the best fitting model 
for spectral ratio studies can only be correlated with 
other stations of the array. A crust with a thickness 
of 35 km is concordant with the ones obtained for 
other stations. In addition, these spectral ratio 
studies favor a multilayered crustal model with rela- 
tively low velocities for P waves. 

Azimuthal deviations for these stations do not 
show any particular pattern (Figure 16) but the arrivals 
from the SE sector are consistently earlier. This 
could be related to a thinning of the low velocity layer 
in the mantle under the plains away from the mountains 
or it could be effectively generated near the source 


region. 


5. TRO. and DEL. ‘stations 


The transverse component of ground motion at 
TRO and DEL showed very large amplitudes; sometimes 
larger than the vertical component (Figure 17). The 
azimuthal deviations show a definite pattern with 
azimuth (Figure 18-a). This pattern of azimuthal 
deviation be modelled with the use of dipping interfaces 
(Niazi, 1966; Zengeni, 1970). The zero crossings of 
azimuthal deviation will determine the direction of 


dipping structure. During the computation of azimuthal 
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Figure V.17 


Three earthquake data recorded by TRO station. 
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Figure V.18 


a) Azimuthal deviations measured as function of 
azimuth, slowness measurements using all 


stations in the array and residuals for TRO. 


b) Azimuthal deviation as a function of distance 
and minimum of ratio of transverse energy to 
the total horizontal energy as a function of 
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azimuth. 
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deviation, clockwise rotation was chosen as a posi- 
tive. Therefore, the second zero crossing will be 

the dip direction. For both DEL and TRO the second 
zero crossing is at about 290° azimuth. A slowness 
measurement using all stations in the array has also 
been shown in Figure 19 (After Gutowski, 1974), in 
which a strong anomaly is clear and it suggests a 
dipping structure under the array. However, this 
anomaly can only be correlated to the azimuthal devia- 
tions observed at TRO and DEL stations. The azimuthal 
deviations showed no definite pattern with distance 
(Figure 18-b). Niazi's tables for determining the dip 
angle are not applicable because the azimuthal devia- 
tions are very large and not regular enough for a simply 
dipping layered model. Since the time residuals at 
these stations are small, any dipping or faulted inter- 
faces must be near the surface and local in nature. A 
refraction study (Richards and Walker, 1959) does not 
indicate any anomaly in crustal structure in a north- 
south profile in central Alberta. If we assume the 
seismometers were positioned and operated with suffi- 
cient accuracy, this anomaly is not a large one because 
it is only observed by two stations and therefore a 
steeply dipping crust to the west under the array cannot 


be postulated. 
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Figure V.19 


Azimuthal deviations measured as a function of 


azimuth and residuals for DEL station. 
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Anomalous results were also noticed on the 
spectral ratios (Figure 20) which have very low am- 
plitudes at these stations. Horizontal components 
on the seismogram (Figure 17) have large amplitudes 
and some phases cannot be correlated with the ones 
on the vertical components. Since horizontal motion 
is generally shear motion, the large arrivals on the 
horizontal components can be related to anomalous P 
to S conversions, possibly in the sedimentary layers. 
Combined spectral ratios for the SE and SW sectors 
have a smooth appearance but they are well separated 
in frequency. The SW events can be modelled by a 
crustal section about 36 km thick but the SE events 
are best modelled by a 10 km thicker crust. Events from 
NW azimuth yielded a spectral ratio with two peaks 
(Figure 20). Splitting of this peak position occurs 
around 290° azimuth. Insufficient data from all 
azimuths makes it difficult to interpret spectral 
ratios in terms of dipping structure. The SE events 
favor either a thicker crust or slower transmission of 
waves whereas the SW events favor a thinner crust or 
faster transmission of waves through the same crustal 
structure. Either of these interpretations is compa- 
tible with the azimuthal deviations. The distance 
between TRO and DEL is about 20 miles and it is clear 


that this anomaly extends over at least this distance; 


pee | : 


gona pnt fit iw Begeloruos ed Fonnes weatie s 

moh sem ‘Egtnos }1or Sons - - Shetccogmnoo mie a 

eft ao elavitxs apzal ‘sag 4 to ktonm soda acter 

7 spoLeamhons of duane - ago e4aenoqmeo tssnosia 
-ereyel yuactnetmibes sd | cb Yihdtewod -ano henewnibg ot 

etodoes we Buse 42a ena trek, nates Ieujoega E : 

badernane Liew ‘sore yearts suid eompuseqas dtoome a 

6 yd belisbom od nso etaeve We ont. sxoneupedt 

adasys #2 ent sud sokds mi VE suoda MOLISE ne sca 

mout ateaeva .sewto retotdd nox OL a Yd’ Besfobom sad a 

| aanaq ows asiw oktes tet3sega & paprany as wy 

BkeoDO wolsLeog Aged aids 2o vadtth La@s (08 owes) 


ae iy : 7 
i 


ne 


Cis. mow? steb tas iokttwenk tidumiss “ORS 


na) 7 
Leasoage sali ot siupitep M4. abit auto 


atneve 88 oT .sxusomete piiqgib 36 amid ak BPR 
to noice tmensay tawole. to seit, seAv ids 6 sodd.ts sari 
10 Savio tettatds & toval atnove we eds esouotiw soveil bs i 


A) fi 
iy 3 
> ae 


BaMoS at eaoltasenqtetnd sagft To teisia .oupsounte ta 


ve 


leseuto empe. sift diguoxss asvaw to noleeimanetd aati 


dériateth GHP .encttsiveh tatduatun was die tds 
wselo ai JE bas eolim Of sede ek gaa brs ont noewsed e 
;sonstelb etd? sasel 36 revo abasoxs ylemons| atdd $eds, 


S enh 
ae int Ag 
re Bi By) 
ligt H 


‘i 
AV OY ah ; 
i y , 4 
y's 5 t 
i 
’ 
is are cia 
3 in | 
( fil nf re He ye sta 
| i 
\ i ee oy 
7 i 
gov). Neen 
te i y 
i ‘ay 
é a iy 
1] Av 
ie in : 
min f , 
ary , , 
L . 
en bee net. BY neo aaa pe 
ary \ GN oy ; ie ait . 
r . hy Waa : 1 s . pa} : ha fe z 
Y y : = i) an iy, | + 
donc it dag ele Shapley ale pear te ela iptina iyelete litre ule aire arena bere ye regal, eee 3 hy Ny 7 
jtagreb A Asm aa ie ae * a ae * wee er Pe we vee” awe | bi: 7 Riis: 
Faerie AVP t) : a ad i.e a wa ea) 
7 : I i | t | ! 
a Feet i Py , 7 i J 
eg a 
? ho i 
: 


saan 


ee EE aes akon its 


Me: Nee 
ae 
hf iy ai! ' aan) ( : i \ A : 7 
: 2 As ge lie CW's 1 
[eat iwte ict = 
ie 7 she f es | aera he eo ‘ i! 
ih * Sie ih 0 ae, ¢* ania ne ay By | fs 
of ae ee * ae ni ie . i 
a a a = t ; 5 as Ail | ite A a arty a o: ’ | 
my fp aie! : ‘ oe . ' wy Nie le a, ; 
wn iis ga. we wrcieh ¥ ae nee | | iS 
; ; " { 
z ne : Baa vile qoieeon rine Wes ‘ ae 4 : | ahd i ; 
Be Pah ii. Nes ’ Gia) 
He Je “" p, ‘ \ ! i ' i 1 , = 
h hall? + Joe f ty ; | 
‘ak jicieiaee | ba | ’ 
ie es Set - r a ; if 
aa ihe is, h aor 
re a i nn ‘ as : 
: Una a ‘ ; 
ome A 15 4h geen A nf vt 
EP ll haan a ae eS : fy 


mavahianar dog’ ern | 
" \ hat Na ae 


Lees b/ : ' é : iy ah j ae eh 
ele nes Wh Giahat en Sa a 4% eh | i) ner 
AN © idl ' tha oC , cay ite mi 
an te F ste , “f ; aa 4 i, ; fn oni } v cane et 9 = tac al He ' "hi f 
a ary " F ia, a 7 : : i aN J mat i , me 
RON ies, eS ame Y ks el oe RN ca beac tain A Rt Lh aie 
Bae Rad gull s. +e an ; Mae) ea i" a nd ea \ ny 1) vor ine i 1, 


Figure V.20 


Long period spectral ratios for TRO and DEL and 


their averages. 
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however, it is not seen at Pine Lake, which is 20 miles 
to the SW. 

From the record character, large horizontal 
motion indicates anomalous conversion of P to S occurs. 
The most probable place will be the low velocity sedi- 
mentary layers. These layers should have highly irre- 
gular structure to account for the conversions. 
Therefore, a steeply dipping or possibly, more compli- 
cated thick sedimentary layers can be postulated to 
explain distortion in the long period waves while the exact 


cause of anomaly cannot be resolved at this time. 


62 Conclusions 


In section III, it has been shown that 35 seconds 
of P coda allows us to obtain sufficient resolution for 
the interpretation of long period spectral ratios. 

Most of the experimental long period spectral ratios 
showed only one dominant peak in the range, 0.05 to 0.2 
Hz but there were some small peaks on each side of the 
Major one. The disappearance of the small peaks in 

some earthquake data may be attributable to the presence 
of additional phases but there were some clear spectra 
from deep events which, even with a larger data window, 
also did not show these small peaks. This can be ex- 


plained by the results obtained in section III, where 
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many model calculations with dipping interfaces showed 
very small or unnoticeable side lobes. Many long 
period experimental spectral ratios presented in this 
section had about 10-30 mHz shift in major peak posi- 
tion from event to event. These variations cannot be 
explained by additional phases alone and may also be 
due to interfaces, dipping as much as 10°. 

By the amount of experimental data taken for 
spectral analysis, we assume an incident plane wave 
which is impinging on the 60-100 km region of the Moho, 
from the station to the opposite direction of wave 
approach. Therefore, assumption of lateral homogeneity 
and no secondary sources is necessary in this part of 
the crust. For the crustal interpretation with a 2-3 km 
sedimentary layer, the model calculations assume that 
about a 20 km region near the station also possesses 
flat or simply dipping layers. Generally, the spectral 
ratios observed at short periods did not match the ones 
obtained by drilling and sonic log data under a particu- 
lar-station’ This would seem to indicate the rapid 
variations and the violation of basic assumptions in this 
method since short period waves average over large 
region of Precambrian crystalline basement. 

EDM long period analyses indicate that the crus- 


tal thickness is about 36 km. The amount of uncertainty 
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is large because there are structural variations with 
azimuth and insufficient amount of uncontaminated data 
for all azimuths. Therefore, there are crustal varia- 
tions in thickness of about 4 km or the velocity 
structure may vary from the mean of 6.5 km/sec. Short 
period spectral ratios for this station are well 
determined as indicated by their duplication on various 
events. The short period interference pattern for this 
station can be approximated sufficiently by a two layer 
sedimentary section. Uncontaminated events show much 
lower amplitudes than this model, but this may well be 
related to the soft layers in the sediments and lower 
velocities in the deep crust. A thin low velocity layer 
in the upper crust cannot be resolved in the frequency 
band available. 

Dependable data for the PIN station shows Similar 
long period spectral ratios to the ones determined for 
the EDM station. RM stations showed completely 
different long period spectral ratios and a model with 
a 35 km crust is suggested. The long period interference 
pattern can be modelled by using a larger number of 
layers than necessary for other stations. TRO and DEL 
data indicates a very heterogeneous structure. Azimuthal 
deviation measurements and record character suggest that 
this anomaly is due to a steeply dipping or more compli- 


cated structure (i.e. with faults) in the upper section. 
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The short period spectral ratios for PIN, TRO and DEL 
show no correlation with each other indicating a 
rapidly Varying structure in the upper crust. Short 
period energy has been concentrated into a small 
frequency band at the RM stations and modelling of 
these spectra is difficult using the assumption of 
plane layering. However, the closely spaced peaks 
indicate a sedimentary section 3-4 km thick, with 
possibly lower velocities in the southern part of the 


array. 
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CHAPTER: £L1 


SYNTHETIC SEISMOGRAMS 


Es Seismograms for P waves 


A time history of the surface motion of a 
layered media responding to a pulse can be obtained 
by convolving the impulse response with a source 
function. Hannon(1964) obtained a time series for 
a number of models using a low frequency pulse. He 
studied the effect of soft layers and variations in 
the angle of incidence on the surface motion. McCamy 
(1967) has made a comparison between experimental and 
theoretical seismograms calculated by the Haskell- 
Thomson matrix method, 

The Haskell-Thomson matrix formulation yields 
a complete plane wave solution for a layered medium. 
Therefore, when a time synthesis is formed, all phases 
will be included provided that the sampling frequency 
is sufficiently small. In this thesis, synthetic 
seismograms were calculated for a number of Alberta 
crustal models by an inverse Fourier transform of the 
frequency components of the ground motion. A fast 
Fourier algorithm was used to obtain the time response 


of a layered medium: 
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is limited by the computational difficulties associated 
with the selection of phases in the partial ray expan- 
Sion. However, a ray theory is more flexible than a 
matrix method, in that the solution is not restricted 

to plane parallel interfaces. A comparison between a 

ray and a wave solution has been carried out to evaluate 
the criteria used in selecting an optimum number of terms 
in the ray expansion (Figure 1). Using a number of test 
parameters outlined in section III, a reasonable match 
between two solutions was obtained for a flat layered 
media. ina ray ‘solution containing 1380 rays for the 
construction of a seismogram the difference in ampli- 
tudes calculated by the two methods was less than 0.1 % 
for the first 3 seconds of time. The accuracy of the 
partial ray expansion decreases with later arrivals 
unless more rays with converted phases are allowed. The 
arrival times of all phases agree and no significant dif- 
ference in amplitudes has been found. A major advantage 
of the partial ray expansion is that all phases are identi- 
fied and may be listed on the Calcomp plot of the synthe- 
tic seismogram. 

A high frequency source is desirable because the 
seismograms will have a high resolution enabling us to 
see differences between the two solutions clearly. 
However, as is often the case with many earthquakes 


the source has a larger period and duration. [In that 
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Figure I.l 


Comparison of seismograms computed by Haskell-Thomson 
and ray theory for a plane P wave impinging on the 
bottom of the layered media (ALTACRT-3) with i= 70° 
(ALPHA = 20°). ‘In: this diagram coding of rays sis 


specified by the equation 
NOFRAY=100000 .JRAY +1000.JCONV + JC 


where JC is an integer which indicates the point of 
incidence at which phase conversions occur for any 
particular ray. JRAY is the ray number from ray 

table (section III) and JCONV is the number of con- 


versions experienced by the JRAY. 
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V- COMPONENT Of PULSE NO.1 FOR ALPHA=20.00° 
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Situation, interference of phases is important and 

in section III, it has been shown that the long 

period spectra obtained either by the ray or wave method 
are in reasonable agreement in position of peaks and 
amplitude. The fast Fourier inversion method of 
obtaining a seismogram, even for a high frequency 

pulse, requires 50 % less computing time when a partial 
ray expansion for the models used. However, the ray 
expansion enables us to interpret the various phases 
easily. For the model shown in Figure 3, there are 
100,000 rays which have 13 segments in the layered 

media but only 0.031 % had amplitudes large enough that 
they had to be included. Rays with a large number of 
segments are very important for the long period spectrum. 
As shown in Figure 3, the number of rays incorporated in 
a dipping model generally is lower. This occurs because 
rays are eliminated by the program when they approach 
the vertex of the wedge. A dipping sedimentary layer 

as much as 10 degrees reduces the number of rays used 

in the solution. The long period spectrum may be un- 
reliable. This restriction is reasonable in a sedimen- 
tary basin because layers seldom have the same dip for 
extended distances and wedging is only common near the 
edge of the basin. A coordinate system has been chosen 


with the +z axis being upward and +x axis in the direction 
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Figure 1.2 


Existing rays and their use for a 4 layered 
crustal model over half space. Dipping model 


indicates only dipping Moho. 
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Figure 1i.3 


Effect of dipping layers on seismograms appearance. 
ALTACRT-1 model was used as in section III. 
Model A is flat layered media 
Model B has a sedimentary section dipping +10° 
Model C has second interface dipping +10° 
Model D has second interface dipping -i0° 


Model E has all layers dipping +10° 
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of phase velocity. The effect of a dipping sedimentary 
layer (Figure 3-B) is greatest when compared to a flat 
layer case (Figure 3-A). A low velocity sedimentary 
section produces large horizontal motion and introduces 
a large phase. Therefore, S motion represents an 
elliptical particle motion with an axis in the quadrant 
formed by negative vertical motion and positive hori- 
zontal motion or positive vertical motion and negative 
horizontal motion. Therefore, elliptical particle 
motion with various angles of tilt in the axis may occur. 
Steeply dipping thicker layers with higher velocity did 
not produce the same effect. Only minor differences 
exist for large time intervals. 

For a flat layered model, a pulse with wavelengths 
larger than layer thickness usually is relatively trans- 
parent to the layer, although minor differences in 
amplitude are noticeable. For example, a number of 
seismograms have been constructed for Model-l and Model-2 
of section V.. For frequencies centered around 1 Hz, 
there is not any significant change but only amplitude 
differences exist. For various models, both the matrix 
method and the ray method have been used simultaneously 
and they produced the same results for wavelengths which 


are 0.015-10 times of the layer thicknesses. 
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1 ee Seismograms for S Waves 


The time response of a layered media to a plane 
S wave incident at the base can also be calculated by 
the Haskell-Thomson matrix method (section II). For 
this case the wave and partial ray expansion techniques 
produced identical results within the limits of compu- 
tational accuracy. In Figure 1, a low frequency 
seismogram (0.5 Hz) calculated by the matrix method 
has been shown for crustal model ALTACRT-3. This model 
has a two layered sedimentary section which is common 
in a continental crust. In this Figure, the main S 
arrival on the vertical component has the second largest 
amplitude but an earlier arrival was identified as a 
conversion at the bottom of the sedimentary layer. 
However, the direct S wave has the largest amplitude 
on the radial component. On vertical component, the time 
difference between converted phase which has largest ampl- 
itude and the pure S phase is delayed about 1 sec and it 
may be as large as 2 seconds depending on thickness and 
velocity of the sedimentary layers, and it was found 
to be not sensitive to dipping interfaces. On the 
vertical component of ground motion, a few other phases 
with significant amplitude exist about six seconds before 
the main phase. These precursors have been identified 


as conversions from S to P at the Moho and at intermediate 
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Figure II.1 


Theoretical seismograms calculated using a 
Haskell-Thomson matrix formulation for a plane 
S wave incident at 30° on the ‘crustal model 
ALTACRT-3. Incident pulse and its spectrum are 


shown in the insert. 


149a 


(SGNO93S) SWIL a 
a ae ee i) pea i 3-0 21-0 o°a, 
WMWLI3dS 3S Wd 3 
R 
2 ' 
Se ee oe a ee ee ee i) 
ae es oe bat th) 88 LHe 
ep) 
2- 
3S Nd 3 
3 WIOW 
a 
= 


“930 O€=1 S LNZOIINI ‘°2H S°O 3S Md WITLYIA 


1SNWS bl 


eg 


——e A a ee en oe 
Eel 
= 


iw 


. 


vj wete Feat 


‘ 
? 


| 2 
Thee bes Ge gwmert, ame caleoud 


(2k COUGS 
2 


# 

3 
> 
& 


Pa 
+ VES " mats 23 f ormilsatin 


he yt 
tT en 


wave smgrie inn t~ab—3u" on the/ae 


il 


NALTACHT T.—tapsiient ib sagen? Ite 


mn 2am the 4 


fees 
. ; a 4 


LIWE 


© 


neEerd. eed 


bs 
» 


¥ 


5 


j 


ll 


sectamen 2) 109s Oe 


ee ee eo 
wa 
ae - 
ed 
4 
7 
a 


£50 


layer. Obviously, these could be used as indicative 
phases for determining crustal thickness. 

The interpretation of S arrival has been always 
a difficult task. Neuman (1930) noticed early longi- 
tudinal arrivals before the main S arrival for an 
earthquake at an epicentral distance of 2500 km. He 
interpreted this as a SP phase and measured 4 seconds 
of time difference between this phase and S5 phase. 
Gutenberg and Richter (1931) could not identify any S 
arrivals for epicentral distances between 9°-15°. 
Byerly (1934) noted the difficulties in identification 
of S phases and he classified S waves as impressive 
(iS) and early phases (eS). Jeffreys and Bullen (1940) 
have noted that S waves are not satisfactorily readable 
up to 25° unless the source is very deep. In Figure 2, 
the amplitude of the vertical and the horizontal com- 
ponent of surface amplitudes has been shown for various 
epicentral distances. The converted phase at the bottom 
of the sedimentary section has the greatest amplitude 
on the vertical component at all distances. On the 
radial component, two phases have very close amplitudes 
for epicentral distances of about 16°. 

The Haskell-Thomson matrix formulation will 
determine the response of a layered media to an input 


pulse for all frequencies when 
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The vertical and horizontal component of 
surface displacement as a function of angle 
of incidence for S and the precursors, SPPP 
and SSSPP from Moho and the base of the 


sedimentary layers. 
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where k is a wave number, w is angular velocity and Vp 
is the maximum velocity in a layered media. This 
condition is equivalent to the phase velocity being 


larger than the maximum velocity in the media. This 


is always true for P waves. That is: 


However, for S waves; 


Ose Ve csc i, 3 (3) 


Therefore, for some epicentral distances conditions 
(1) and (2) are violated. The maximum wave number 

to be evaluated falls beyond the poles of matrix for- 
mulation in the complex k plane (pene vi a Wes) and 
they are imaginary for P waves, producing a damped 
motion in frequency. In the matrix elements, terms 
containing a sin or cos will.be sinh and cosh. 
Therefore, for these cases, the response can be 
calculated for wave numbers that occur before the 
first pole, k = w/Vp which corresponds to C = Vj. 


The case where C < V., produces a phase dif- 


P 
ference of +90° between the vertical and the horizontal 
components of S motion which corresponds to an ellip- 


tical particle motion, the axis of the ellipse being 
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vertical and radial. This arises because the numerator 
of the transfer function for the horizontal component 
will be imaginary and the numerator of the transfer 
function for the vertical component will.be real. The 
transfer function for S waves oscillates rapidly and 
produces a complicated particle motion. As far as we 
know, there have been no studies on the change of fre- 
quency of S waves with epicentral distance. But Nuttli and 
Whitmore (1962) have investigated the particle motion 
of S waves. They found that the interpretation can be 
done best when the phase velocity is larger than the 
maximum velocity in the crust. 


The condition where C < V, in a layer corres- 


2) 
ponds to the propagation of "trapped" modes for P 

waves but for the solution to exist the phase velocity 
must be greater than the S wave velocity a any layer. 
The ray theory approach should produce the same 

results as the matrix method. However, the ray approach 
used here does not predict the frequency characteristics 
of the surface motion. In Figure 3, seismograms still 
retain their high frequency character (10 Hz). For.an 
angle of incidence of 70°, ‘the phase velocity is less 
than the P wave velocities in a few layers. However, 

in Figure 3, the seismograms have a flat region 


between 7-11 seconds. These regions are due to rays 


which do not reach the surface since C < Vp in the 
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Figure 11.3 


Theoretical seismogram for a plane S wave 
incident at an angle of 70°. . Among 10976 


created rays, 1613 of them were used. 
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model. In the ray method returning rays have not 
interacted with upgoing rays. However, this has 
been done in the wave solution producing an attenu- 
ation of S waves with a converted segments. 
Precursors to the S or shear phase due to the 
presence of a crustal layer can introduce serious 
errors in interpretation if only the vertical com- 
ponent is available. The error may be as much as 
1-6 seconds at all distances and is of particular 
-t.,, times are used to determine 


Sy 


epicentral distances, first motion studies and in the 


importance when t 


study of premonitory phenomena using the ratio t./tp 


for earthquake prediction. 
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